Dw = 04188

RR312

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

ARCHIVES

ÉTUDE, RÉSOLUTION
ET OPTIMISATION
DE MODÈLES
MACROÉCONOMIQUES

Daniel GABAY
Pierre NEPOMIASTCHY
M'hamed RACHDI
Alain RAVELLI

Rapport de Recherche N° 312

Juin 1978

68p
# ÉTUDE, RÉSOLUTION ET OPTIMISATION DE MODELES MACROÉCONOMIQUES

Daniel Gabay, Pierre Nepomiastchy, M'hamed Rachdi, Alain Revelli
IRIA/LABORIA

## Résumé :

Après avoir décrit la structure générale d'un modèle macroéconomique dynamique, non linéaire et déterministe, on montre comment on peut renuméroter les équations du modèle pour rendre sa jacobienne quasitriangulaire. Des algorithmes de résolution et d'optimisation adaptés aux modèles quasitriangulaires sont proposés. Les comparaisons numériques avec les méthodes classiques sont faites sur un modèle réel de l'INSEE.

## Abstract :

This paper describes the general structure of a dynamic, non linear and deterministic macroeconomic model, then shows how a renumbering of the model equations can put the jacobian matrix in a quasitriangular form. Adapted to quasitriangular model algorithms for solving and optimizing the model are proposed. Numerical comparisons with classical methods are performed using a real model of INSEE.
ETUDE, RESOLUTION ET OPTIMISATION

DE MODELES MACROECONOMIQUES

D. GABAY, P. NEPOMIASTCHY, M. RACHDI, A. RAVELLI

IRIA - LABORIA
B.P. 105
78150 Rocquencourt

AVANT-PROPOS.

Le travail exposé dans cet article est le résultat de recherches
sur les modèles macroéconomiques menées au LABORIA depuis deux ans
environ. En juillet 1977, les chercheurs concernés ont été regroupés
au sein d'un nouveau projet LABORIA baptisé *Résolution et optimisation
de modèles macroéconomiques*.

Par ailleurs, en avril 1977, le LABORIA a décidé de lancer une
action pour la construction d'un logiciel spécialisé pour la construc-
tion et l'utilisation d'un modèle macroéconomique. Ce logiciel, basé
sur la programmation modulaire, sera chargé d'effectuer et de coordonner
des tâches aussi diverses que la gestion d'une banque de données écono-
miques, l'estimation des paramètres du modèle, la simulation etc...
Vu le caractère fondamentalement pluridisciplinaire des travaux impliqués,
la construction de ce logiciel sera effectuée par des économistes,
des statisticiens, des mathématiciens et des informaticiens de plusieurs
organismes réunis au sein du *Club MODULECO* animé par le LABORIA (voir
[1]).

En plus de leur intérêt scientifique propre, les recherches décrites
dans ce rapport serviront de base à la programmation des modules de
simulation et d'optimisation du logiciel MODULECO.
# I. INTRODUCTION.

Un modèle macroéconomique est un système d'équations aux différences dont la solution est supposée décrire le comportement d'une économie nationale. Ces modèles sont de plus en plus utilisés en France pour l'élaboration de prévisions économiques (préparation du Budget) et pour la génération de scénarios du Plan.

On constate ces dernières années une augmentation de la taille de ces modèles. Par exemple, le modèle DMS [2], l'un des modèles les plus récents construits par l'INSEE, comporte plus de 1000 équations pour chaque période. Dans ces conditions, on conçoit aisément l'importance que peut présenter la mise au point de méthodes efficaces pour la résolution et l'optimisation de ces modèles.

Jusqu'à une date très récente, l'absence de recherches pluridisciplinaires dans le domaine (c.a.d. rassemblant des économistes et des mathématiciens) faisait que ces problèmes étaient résolus par utilisation brutale de packages généraux de recherche opérationnelle. Ces codes ne sont évidemment pas adaptés à la spécificité des modèles macroéconomiques ; on montre dans ce papier que l'adaptation des algorithmes généraux à cette structure spécifique peut diviser le temps de calcul par un facteur 10 ou plus. Ce fait est aujourd'hui bien perçu et au dernier congrès du National Bureau of Economic Research, une session entière a été consacrée à la résolution et à l'optimisation de modèles macroéconomiques. A cette session, trois papiers, Drud [3], Lasdon [4], Népomiastchy-Ravelli [5], traitaient de l'adaptation de méthodes connues (méthodes de Newton, de l'état adjoint, du gradient réduit) à la structure particulière qu'ont les matrices de Jacobi du système d'équations décrivant les modèles (matrices très creuses et quasi-triangulaires inférieures).

Nous donnons ici un panorama aussi général que possible des techniques de résolution (chapitre III) et d'optimisation (chapitre IV) de modèles macroéconomiques (utilisées ou non dans des logiciels existants). Nous parlerons entre autres des méthodes de résolution et d'optimisation développées au LABORIA. Ces méthodes exploitent la structure quasi-triangulaire du jacobien du modèle. Celui-ci nous a parfois été donné par ses auteurs déjà sous une forme quasi-triangulaire (modèle PIMPON de la Direction de la Prévision [6] et modèle STAR de l'INSEE [7]). Ce n'est
pas toujours le cas, et à l'occasion de l'étude du modèle de R.C. Fair [18],
nous avons mis au point un algorithme de renumérotation automatique des
équations et des variables afin de rendre le modèle le plus triangulaire
possible. Cet algorithme est décrit au chapitre II.

Le traitement des matrices creuses pour la méthode de l'état-adjoint
ainsi que l'objectif final de cette étude (implémentation de nos algorithmes
dans la bibliothèque MODULECO) nous ont conduit à faire des recherches en
programmation, qui sont décrites au chapitre V, et dont l'objet est de
rendre nos programmes aussi efficaces, généraux et faciles à utiliser que
possible.

Une autre utilisation importante des modèles macroéconomiques est l'étude de la sensibilité de l'économie à une variation des variables d'entrée (prix du pétrole par exemple) ou des coefficients du modèle, ce qui revient à calculer les **multiplicateurs dynamiques** du modèle. Ces derniers sont en général calculés par une série de simulations, comme pour les multiplicateurs du modèle DMS [8]. Une étude plus systématique et surtout plus globale peut être faite par une linéarisation du modèle autour d'une trajectoire de référence, puis du calcul des valeurs propres de la matrice du transfert. C'est l'approche de Oudet pour l'étude de la dynamique du modèle STAR [9]. Cette approche nous paraît très intéressante mais, ayant été faite à l'aide de packages généraux d'analyse spectrale, est coûteuse pour les gros modèles. Il semble qu'il y ait là un besoin : trouver des méthodes d'analyse spectrale adaptées à la structure particulière des modèles macroéconomiques. Pour une étude générale de l'analyse des modèles, nous renvoyons le lecteur à Deleau-Malgrange [10].

Les mécanismes fondamentaux de l'économie étant difficiles à saisir
à l'aide d'un modèle de 1000 équations ou plus, les économistes sont très
intéressés par la construction de "modèles réduits" appelés maquettes,
comme la maquette MINIDMS [11] du modèle DMS. Ces maquettes sont souvent
construites par une détermination a priori des variables importantes du
gros modèle. Une étude plus systématique, celle de Dufour-Gilles [12]
revient à mesurer le degré d'interdépendance de chaque couple de variables
à l'aide d'une norme tirée de la théorie de l'information et d'éliminer
toutes les connexions au-dessous d'un certain seuil.
Enfin, un domaine que nous n'aborderons pas ici est celui des modèles stochastiques. Nous renvoyons le lecteur à Bensoussan-Hurst-Naslund [13] pour la théorie générale des systèmes stochastiques, à Athans [14] pour la validité des modèles en fonction du degré d'incertitude, à Kendrick [15] et Wall-Wescott [16] pour l'application de la théorie du contrôle stochastique en macroéconomie et, enfin, à l'étude bibliographique de Pau [17] sur les modèles macroéconomiques stochastiques.
## II. ETUDE DE LA STRUCTURE DU MODELE.

### II.1. Notations.

On considère un modèle dynamique en temps discret pour des périodes $t$ allant de 1 à $T$. Dans les modèles courants, le nombre $T$ varie de 5 à 30 et la période peut être un trimestre ou un an.

Les variables déterminées par le modèle, appelées variables d'état par les automaticiens et les utilisateurs de la théorie du contrôle optimal, sont appelées **variables endogènes** par les économistes et nous adopterons ici cette terminologie. Soit alors $n$ le nombre de ces variables endogènes; on note $x_t^i$ la valeur de la variable $i$ à la période $t$, $x_t$ le vecteur $\{x_t^1, \dots, x_t^n\}$, $x_t \in \mathbb{R}^n$, et $x$ le vecteur $\{x_1, \dots, x_T\}$ $x \in \mathbb{R}^N$ avec $N = nT$.

Soient $r$ le nombre de **variables de contrôle** (ou variables de décision) du modèle, $u_t^i$ la valeur du contrôle $i$ à la période $t$, $u_t$ le vecteur $\{u_t^1, \dots, u_t^R\}$, $u_t \in \mathbb{R}^R$ et $u$ le vecteur $\{u_1, \dots, u_T\}$, $u \in \mathbb{R}^T$ avec $R = rT$.

Soient $p$ et $q$ les retards maxima apparaissant dans le modèle sur (respectivement) les variables endogènes et les variables de contrôle.

Avec ces notations, le modèle est entièrement décrit par le système d'équations :

$$ (1a) \qquad \begin{gathered} k_i x_t^i = f_t^i(x_t, x_{t-1}, \dots, x_{t-p}, u_t, u_{t-1}, \dots, u_{t-q}), \\ i = 1, \dots, n \quad \text{et} \quad t = 1, \dots, T \end{gathered} $$

$$ (1b) \qquad x_t^i \text{ donné pour } i = 1, \dots, n \text{ et } t = 0, -1, -2, \dots, -p+1 $$

$$ (1c) \qquad u_t^i \text{ donné pour } i = 1, \dots, r \text{ et } t = 0, -1, \dots, -q+1 $$

Dans la notation (1), les variables économiques qui sont considérées comme des données du modèle, c'est à dire les **variables exogènes**, sont incluses dans les fonctions $f_t^i$.

Les $k_i$ apparaissant dans (1a) sont des paramètres donnés valant 0 ou 1.

Dans un modèle macroéconomique, les équations dites **équations de comportement** sont estimées comme **relations liant une variable**
(variable expliquée) à plusieurs autres (variables explicatives) ;
l'écriture $x_t^i = f_t^i(x,u)$ est dans ce cas naturelle et correspond au choix
$k_i = 1$ dans (1a). D'autre part, parmi les autres équations du modèle
(équations de définition et équations comptables), beaucoup peuvent être
résolues analytiquement par rapport à la variable portant le même
numéro que l'équation (c'est en général une simple question de numérota-
tion), ce qui correspond encore à $k_i = 1$. Cette écriture $x_t^i = f_t^i(x,u)$
est plus commode et plus économique pour plusieurs algorithmes de
résolution et d'optimisation (cf. chapitres III et IV). Comme elle n'est
pas toujours possible (¹), nous avons adopté la notation générale (1a).

Par conséquent $k_i = 1$ si l'équation i est résolue en $x_t^i$ et $k_i = 0$ ou bien si $x_t^i$ n'intervient pas dans l'équation i ou bien si $x_t^{i}$ intervient mais si l'équation i n'est pas résoluble analytiquement en $x_t^i$.

Nous supposons connu un sous-ensemble D de $\mathbb{R}^R$, domaine des valeurs "raisonnables" du contrôle u; nous supposons que les fonctions $f_t^i$ sont deux fois continûment différentiables par rapport à toutes leurs variables et que, au moins pour tout $u \in D$, le système (1) a une solution et une seule, notée $x(u)$. Le modèle (1) est alors équivalent à la donnée d'une application $u \to x(u)$ de D dans $\mathbb{R}^N$.

Nous noterons K la matrice diagonale dont le ième élément diagonal est $k_i$ :

$$
(2) \quad K = \operatorname{diag}(k_i)
$$

Nous noterons $d_t$ le vecteur de $\mathbb{R}^m$ (avec $m = pn + (q+1)r$) :

$$
(3) \quad d_t = \{x_{t-1}, x_{t-2}, \dots, x_{t-p}, u_t, u_{t-1}, \dots, u_{t-q}\}
$$

qui lors de la résolution du modèle pour la période t, contient toute
l'information connue.

Avec les notations (2) et (3), on voit que résoudre le modèle pour
la période t, c'est résoudre le système algébrique de dimension n :

$$
(4) \quad K x_t = f_t(x_t, d_t)
$$

où K, $f_t$ et $d_t$ sont connus.

<footnote>(1) pour le modèle STAR (voir fig.2, page 13), 136 de ses 139 équations s'écrivent sous la forme $x_t^i=f_t^i(x,u)$ ; pour le modèle de FAIR (voir fig.3, page 14), 82 de ses 83 équations sont de cette forme.</footnote>
II.2. <u>Notion de variable de bouclage.</u>

Soit I un sous-ensemble de l'ensemble d'indices {1,2,...,n}.
On appelle S(I) le système algébrique de dimension m = card (I) obtenu
en ne conservant dans le système (4) que les équations i pour i ∈ I et
en considérant que les variables x_t^j pour tout j ∉ I sont données. Bien
entendu, si I = {1,2,...,n} , S(I) coïncide avec (4).

Pour tout I, on appelle _matrice d'incidence_ du système S(I) de
dimension m la matrice carrée de même dimension dont l'élément e_{ij} est
égal à 1 si la variable x_j^t apparaît dans l'équation i pour au moins l'une
des périodes t (¹) et e_{ij} est égal à 0 dans le cas contraire.

On dit qu'un système S(I) est _récursif_ si sa matrice d'incidence
vérifie :

$$
\begin{align*}
e_{ii} &= 1, \quad \forall i \in I \\
e_{ij} &= 0, \quad \text{pour tout } j > i
\end{align*}
$$

et si, pour tout i ∈ I, l'équation i est résolue par rapport à x_t^i,
c'est à dire si k_i = 1.

Il est clair que la résolution d'un système récursif se ramène
simplement à l'application successive des formules définissant les
équations. On notera que le système global (4) est récursif si et seule-
ment si il est _triangulaire inférieur_ avec k_i = 1 pour tout i.

On dit qu'une variable *j* est une *variable de bouclage* du système (4) ou bien si $k_j = 0$ ou bien si la $j$-ème colonne de la matrice d'incidence de (4) possède au moins un élément non nul au-dessus de la diagonale principale, c'est à dire que *j* est une variable de bouclage si et seulement si :

$$
(6) \quad k_j \prod_{i=1}^{j-1} (1 - e_{ij}) = 0
$$

Une équation portant le même numéro qu'une variable de bouclage est
appelée _équation de bouclage_.

<footnote>(1) plus précisément, la variable $x_t^j$ apparaît dans l'équation i à t si et
seulement si la dérivée par rapport à $x_t^j$ de l'expression
$k_i x_t^i - f_t^i(x_t, d_t)$ n'est pas identiquement nulle.</footnote>
Nous appellerons s le nombre de variables de bouclage, B l'ensemble des indices des variables de bouclage et $\bar{B}$ l'ensemble des indices des variables qui ne sont pas des variables de bouclage.

On notera que les définitions ci-dessus impliquent directement
le résultat suivant, fondamental pour la suite :

**Théorème 1.** : le système $S(\bar{B})$ est récursif.

En d'autres termes, pour des valeurs connues des variables de
bouclage, le calcul des autres variables endogènes est une simple
application des formules définissant le modèle. Il est clair, et on le
verra en détails au chapitre III, que le système (4), qui est de dimension
n, peut se ramener à un système de dimension s dont les inconnues sont
les variables de bouclage et dont l'évaluation des équations revient à
la résolution d'un système récursif de dimension n-s.

On voit intuitivement, et nous en aurons confirmation en III, que plus s est petit et plus le modèle est facile à résoudre. Par ailleurs, il est bien évident que la matrice d'incidence dépend de la numérotation des variables et des équations, et il en est donc de même du nombre de variables de bouclage. C'est la raison pour laquelle nous avons mis au point un algorithme, décrit en II.4, de renumérotation des variables et des équations pour diminuer le plus possible le nombre des variables de bouclage. Cet algorithme, appliqué au modèle de Fair [18] a permis de diviser par cinq le nombre de variables de bouclage.

Si, éventuellement après renumérotation, s est petit devant n, la
matrice d'incidence du modèle est quasi-triangulaire inférieure, c.a.d.
que la plupart de ses éléments non nuls se trouvent au-dessous de la
diagonale principale et on dit alors que le modèle est quasi-triangulaire ;
la plupart des modèles macroéconomiques peuvent être décrits sous une
forme quasi-triangulaire. Nous donnons ci-dessous le degré de récursivité
des modèles que nous avons analysés :
<table><thead><tr><th>Modèle</th><th>Nombre d'équations (n)</th><th>Nombre de variables de bouclage (s)</th></tr></thead><tbody><tr><td>ANDOMINI [19]</td><td>4</td><td>1</td></tr><tr><td>PIMPON [6]</td><td>14</td><td>2</td></tr><tr><td>FAIR [18]</td><td>83</td><td>6</td></tr><tr><td>STAR [7]</td><td>139</td><td>3</td></tr><tr><td>MPS *</td><td>400</td><td>20</td></tr><tr><td>DMS [2]</td><td>plus de 1 000</td><td>moins de 100</td></tr></tbody></table>

### - Table 1 -

Les variables de bouclage ont été obtenues "à la main" pour le modèle ANDOMINI ; elles ont été obtenues par l'algorithme décrit en II.4 pour le modèle de Fair, qui contenait au départ 30 variables de bouclage ; elles ont été fournies par les auteurs du modèle pour PIMPON et STAR et aucune renumérotation n'a été faite pour le modèle DMS.

### II.3. Différentes parties d'un modèle.

Soit donc un modèle qui, éventuellement après renumérotation des équations et des variables, contient s variables de bouclage. L'application de l'algorithme décrit au paragraphe II.4 permet de "découper" le modèle en quatre blocs : le prologue, le cœur, les équations de bouclage et l'épilogue ; ce "découpage" facilite la résolution du modèle. Pour certains modèles, le prologue ou l'épilogue peuvent être vides ; ce ne sont que des cas particuliers et toutes les méthodes décrites ici restent valables (avec $l = 0$ ou $m+s = n$).

Le modèle commence par un nombre $l$ d'équations ne contenant pas les variables de bouclage. Par définition de ces dernières, le système formé par ces $l$ équations, appelé prologue du modèle, est totalement récursif :

$$ (7) \quad \underline{\text{Prologue}} : \left\{ \begin{array}{l} x_t^1 = f_t^1(d_t) \\ x_t^i = f_t^i(x_t^1, \dots, x_t^{i-1}, d_t) \\ x_t^l = f_t^l(x_t^1, \dots, x_t^{l-1}, d_t) \end{array} \right. $$

<footnote>* Ces résultats concernant le modèle MPS de la Federal Reserve Bank ont été obtenus au moment de la correction de l'article, en utilisant l'algorithme décrit en II.4.</footnote>
Ce système peut être résolu directement, indépendamment des autres équations du modèle. Pour donner un exemple d'équation du prologue, dans le modèle STAR la population a été estimée non en valeur, mais par le taux de croissance démographique. Une équation du prologue de STAR est donc $POP_t = \alpha POP_{t-1}$. C'est évidemment le cas le plus simple.

Après le prologue, nous trouvons dans le modèle les "vraies" variables endogènes $x_t^{\ell+1}, \dots, x_t^m$ et les variables de bouclage $x_t^{m+1}, \dots, x_t^{m+s}$. Nous appellerons $y_t$ le vecteur $\{x_t^{m+1}, \dots, x_t^{m+s}\}$ de ces variables de bouclage. Par définition de ces dernières, le système d'équations $(\ell+1)$ à $(m)$ est récursif pour des valeurs données des variables de bouclage, le prologue ayant été résolu au préalable. Cet ensemble de $m-\ell$ équations constitue le **coeur** du modèle et peut s'écrire :

$$ (8) \quad \text{\underline{Coeur}} : \quad x_t^i = f_t^i(x_t^1, \dots, x_t^{i-1}, y_t, d_t), \quad i = \ell+1, \dots, m $$

Les équations suivantes sont les équations de bouclage qui, par définition, ne sont pas récursives et, en particulier, ne sont pas nécessairement résolues par rapport à $x_t^i$. Elles s'écrivent :

$$ (9) \quad \text{\underline{Equations de bouclage}} : \quad k_i x_t^i = f_t^i(x_t^1, \dots, x_t^{m+s}, d_t), \quad i=m+1, \dots, m+s $$

Une fois que le système d'équations formé par (7), (8), (9) est résolu, un certain nombre de **variables de sortie** se déduisent directement de la connaissance des variables précédentes par des équations, généralement très simples, qui forment l'**épilogue** du modèle :

$$ (10) \quad \text{\underline{Epilogue}} : \quad x_t^i = f_t^i(x_t^1, \dots, x_t^{i-1}, d_t), \quad i = m+s+1, \dots, n $$

Par exemple, le coeur de STAR détermine les importations et les exportations et une équation de l'épilogue est $BACO_t = EXP_t - IMP_t$, la balance commerciale intervenant dans le critère.

Pour les modèles STAR, FAIR et MPS les dimensions respectives des sous-systèmes (7), (8), (9) et (10) sont les suivantes :
<table><caption>- Table 2 -</caption><thead><tr><th></th><th>FAIR</th><th>STAR</th><th>MPS</th></tr></thead><tbody><tr><td>Prologue</td><td>2</td><td>58</td><td>184</td></tr><tr><td>Coeur</td><td><b>68</b></td><td>73</td><td>162</td></tr><tr><td>Equations de bouclage</td><td>6</td><td>3</td><td>20</td></tr><tr><td>Epilogue</td><td>7</td><td>5</td><td>34</td></tr><tr><td>Total</td><td>83</td><td>139</td><td>400</td></tr></tbody></table>

Dans le cas particulier du modèle STAR, les 48 premières équations du prologue ne contiennent pas le contrôle et, de plus, ne font intervenir avec des retards que les 48 premières variables endogènes. Ces 48 équations forment donc un bloc qui peut être résolu une fois pour toutes juste après la lecture des données. Ces variables, appelées *variables prédéterminées* une fois calculées, on pourra ensuite travailler sur un modèle de 91 équations seulement.

La matrice d'incidence du modèle, définie au début du paragraphe précédent, a donc, par définition des quatre parties du modèle la configuration suivante :

<figure><img src="image_1.png" /><figcaption>- Figure 1 -</figcaption></figure>
dans laquelle la partie non hachurée de la matrice est celle qui ne contient que des zéros.

Les figures 2 et 3 donnent, à titre illustratif, les matrices d'incidence des modèles STAR et FAIR ; dans les deux cas, les zéros de la matrice ont été remplacés par des blancs et de plus, les 1 de la matrice de STAR ont été remplacés par des points.

Il faut ajouter que certains modèles de grande taille, ainsi d'ailleurs que tous les modèles internationaux interrégionaux, ont de plus une structure par blocs, chaque bloc ayant la configuration donnée en figure 1.

Nous utiliserons au chapitre IV la matrice de Jacobi $F_t = ((\partial f_t^i / \partial x_t^j))$ du système (4). On vérifiera que le nombre d'éléments non identiquement nuls de la matrice $F_t$ est égal au nombre total de 1 de la matrice d'incidence diminué du nombre de 1 de la diagonale principale de celle-ci qui sont situés à un rang correspondant à un coefficient $k_i = 1$. Nous donnons ci-dessous des statistiques sur les matrices de Jacobi des modèles STAR, FAIR et MPS, étant entendu que le modèle STAR est ici le modèle de 91 équations obtenu en enlevant le bloc des 48 premières équations mentionné ci-dessus :

<table><thead><tr><th></th><th>STAR</th><th>FAIR</th><th>MPS</th></tr></thead><tbody><tr><td>Dimension de la matrice F<sub>t</sub></td><td>91</td><td>83</td><td>400</td></tr><tr><td>Nombre d'éléments non identiquement nuls de F<sub>t</sub></td><td>173</td><td>218</td><td>865</td></tr><tr><td>Pourcentage d'éléments non identiquement nuls de F<sub>t</sub></td><td>2,09 %</td><td>3,16 %</td><td>0,54 %</td></tr></tbody></table>

- Table 3 -

Ces chiffres montrent de façon éloquente à quel point les matrices de Jacobi des modèles STAR, FAIR et MPS sont creuses et justifient les techniques de traitement des matrices creuses exposées au chapitre V.
<figure><img src="image_2.png" /><figcaption>- Figure 2 -</figcaption></figure>

Matrice d'incidence du modèle STAR.
<figure><img src="image_3.png" /><figcaption>- Figure 3 -</figcaption></figure>

Matrice d'incidence du modèle FAIR.
**Remarque** : il ne faut pas déduire de la table 3 que le modèle STAR ne compte en moyenne que $173/91 = 1,91$ variables extradiagonales par équations, car nous n'avons comptabilisé ici que les variables endogènes non retardées. En ajoutant les variables endogènes retardées, les variables de contrôle et les variables de la diagonale, on aboutit pour le modèle STAR à une moyenne de $336/91 = 3,69$ variables par équation.

## II.4. Renumérotation des équations du modèle.

Il est clair qu'en général, la matrice d'incidence obtenue à partir d'un modèle donné brut, n'aura pas du tout la structure de la figure 1 qui dépend de façon essentielle de l'ordre des équations. Il suffit en effet, en général, d'intervertir deux équations pour faire apparaître un nouvel élément non nul au-dessus de la diagonale principale.

Des études ont déjà été faites pour décomposer l'ensemble des variables en blocs d'interdépendance ce qui se traduit sur la matrice d'incidence par l'apparition de m blocs diagonaux du type de la figure 4 :

<figure><img src="image_4.png" /><figcaption>- Figure 4 -</figcaption></figure>
dans laquelle la partie non hachurée de la matrice est encore celle qui
ne contient que des zéros.

D'autre part, il existe des algorithmes pour déterminer de tels
blocs diagonaux de taille minimale (voir Drud [3] ou Hellerman et
Rarick [20], et qui peuvent encore s'obtenir par élimination ou "tearing"
[voir Steward [21]]). La résolution se fait alors récursivement pour
chaque bloc d'interdépendance par une méthode itérative globale à
l'intérieur de ces blocs.

On aborde dans ce chapitre le problème de la détermination
systématique des variables de bouclage. Ceci signifie une "réduction"
de la structure d'autant plus que le sous-ensemble des variables du
modèle dépendant des variables de bouclage (systèmes (8) et (9)) constitue
souvent à lui seul un bloc d'interdépendance irréductible (c'est le cas
pour le modèle de FAIR).

Les exemples de modèles donnés plus haut montrent l'intérêt d'une
telle étude.

D'autre part, la méthode de détermination des variables de bouclage
proposée peut être adaptée à chaque bloc d'interdépendance. Elle
semble permettre ainsi une bonne compréhension du modèle. Pour l'utilisa-
tion de la méthode de Newton par exemple, (décrite au chapitre III), on
obtient une simplification par rapport à une méthode de Newton globale,
grâce à la réduction considérable de la taille du système par une simple
opération de réordonnancement des équations (opération qui peut être
effectuée une fois pour toutes et avant toute simulation).

Enfin, on peut remarquer qu'il ne s'agit là que d'une méthode
qualitative d'analyse d'interdépendance entre les variables et non pas
d'une méthode quantitative qui tiendrait compte du poids de chaque
variable dans chaque équation avec réduction par niveau d'influence
globale (voir Dufour et Gilles [12]).

L'algorithme est basé sur l'analyse d'interdépendance entre chaque
variable directement sur la matrice d'incidence introduite au paragraphe
II.2.
L'analyse peut être rendue plus claire par l'étude du graphe correspondant au système à résoudre qui est obtenu en associant deux variables $x_i$ et $x_j$ par un chemin orienté (ou fléché) de $x_i$ à $x_j$ si la variable $x_i$ intervient dans l'équation j qui définit la variable $x_j$. A titre d'exemple, on a reproduit à la figure 8 le graphe associé au sous-modèle conjectural-structural [22] dont l'étude faite à la main sert plus loin d'illustration pour notre algorithme. La notion de variable de bouclage est alors étroitement liée à la notion de *boucle* à l'intérieur de ce graphe. Le système ne peut en effet être résolu globalement récursivement que si le graphe associé ne contient pas de boucle. Les variables de bouclage peuvent être alors interprétées comme des singularités dans cette récursivité. Mais il est clair que la mise en place d'un algorithme directement lié à l'étude du graphe n'est pas aisée, ce graphe devenant rapidement indéchiffrable lorsque la taille du modèle grandit.

La matrice d'incidence est beaucoup plus facile à déterminer et peut être calculée soit à la main, soit de manière automatique par une méthode analogue à celle de la détermination de la matrice jacobienne introduite au paragraphe II.3 (différences finies ou programmation par analyse syntaxique des instructions par exemple) puisque, mis à part le remplissage adéquat de la diagonale principale, les éléments identiquement égaux à 1 de la matrice d'incidence correspondent aux éléments non identiquement nuls du jacobien.

Il est nécessaire d'abord de préciser l'importance du tableau K introduit au paragraphe II.1 qui d'entrée indique pour chaque équation si celle-ci est résoluble analytiquement par rapport à la variable qu'elle définit. Ceci est essentiel (et on le précisera au chapitre III) pour la résolution par une méthode itérative. Cherchant d'autre part à minimiser le nombre des variables de bouclage, on peut a priori chercher à minimiser le nombre des coefficients $k_i$ nuls pour lesquels, par définition, la variable $x_t^i$ est une variable de bouclage.

Le tableau K, comme on va le voir, peut différer sensiblement, suivant que l'on décide de faire une analyse grossière ou fine de chaque équation du système à résoudre.
Ainsi, l'équation $f(x_1,x_4,x_8) = 0$ peut en général s'écrire
$x_1 = f_1(x_4,x_8)$ ou $x_4 = f_4(x_1,x_8)$ ou $x_8 = f_8(x_1,x_4)$. En choisissant
l'une des trois variables $x_1, x_4$ ou $x_8$ comme la variable définie
par l'équation $f(x_1,x_4,x_8) = 0$, on peut donc, si l'une des trois
expressions $f_1, f_4$ ou $f_8$ est connue analytiquement, réordonner les
équations de façon à ce que la variable $x_1, x_4$ ou $x_8$ suivant le cas soit
définie par cette expression et ainsi faire apparaître un nouveau coeffi-
cient 1 dans le tableau K. D'autre part, si aucune des expressions $f_1,$
$f_4$ et $f_8$ n'est connue analytiquement, on peut imaginer (dans un programme
de résolution récursive ou par exemple les variables $x_1$ et $x_4$ ont déjà
été calculées) d'introduire un programme de résolution à une variable de
l'équation $f(x_1,x_4,x_8) = 0$ qui permettrait donc de calculer la valeur de
$x_8 = f_8(x_1,x_4)$ ce qui revient encore d'une manière détournée à introduire
un nouveau coefficient 1 dans le tableau K.

On a un problème analogue quand un certain nombre de variables
du modèle sont solutions (à l'intérieur d'un bloc) d'un système linéaire
qui n'est pas résolu analytiquement dans le modèle. Ainsi dans le modèle
de FAIR [18] , après avoir déterminé les variables CG<sub>t</sub>, DDF<sub>t</sub>, BORR<sub>t</sub>,
SAVB<sub>t</sub>, LF<sub>t</sub> et SECR<sub>t</sub>, les cinq équations 45, 61, 62, 64 et 70 forment un
système linéaire de cinq équations à cinq inconnues : BR<sub>t</sub>, A<sub>t</sub>, DDB<sub>t</sub>,
LBVBB<sub>t</sub>, DDH<sub>t</sub>.

Ces équations numérotées de 1 à 5, pour les variables appelées $x_1, x_2, x_3, x_4$ et $x_5$, s'écrivent :

1. $x_1 = ax_3$

2. $x_2 = -x_5 + b$

3. $x_3 = x_5 + c$

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

5. $0 = -x_4 - x_2 + e$

Dans une analyse grossière, équation par équation du modèle, l'équa-
tion 5 (équation 70 du modèle), qui dans FAIR est la seule à ne pas
définir explicitement une des variables endogènes, peut donc a priori
et de façon naturelle, déterminer la variable restante : RBILL<sub>t</sub>, mais
avec l'introduction d'un coefficient $k_{70} = 0$ (RBILL<sub>t</sub> est alors automa-
tiquement une variable de bouclage).
D'autre part, la résolution du système 1, 2, 3, 4 et 5 permet
d'écrire :

1' $ax_3 = b + c + d - e$

2' $x_1 = ax_3$

3' $x_2 = -x_3 + b + c$

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

5' $x_5 = x_3 - c$

Ainsi, l'équation 5 (équation 70), devenue 5', a été "utilisée"
pour déterminer la variable $DDH_t$ et le coefficient $k_{70}$ devient égal
à 1. Parallèlement, l'équation 8 dans le modèle de FAIR qui précédemment
déterminait de façon naturelle la variable DDH et qui s'écrit :

$$
\log(\mathrm{DDH}_t) = -0.0255 * \log(\mathrm{RBILL}_t) + 0.432 * \log(\mathrm{YH}_t) + f
$$

peut aussi s'écrire :

$$
0.0255 * \log(\text{RBILL}_t) = - \log(\text{DDH}_t) + 0.432 * \log(\text{YH}_t) + f
$$

et ainsi servir à déterminer la variable RBILLt, le coefficient k8
restant égal à 1.

C'est de toute façon ensuite notre algorithme qui permet de déter-
miner si RBILL_t est une "vraie" variable de bouclage ou non.

Il est clair que ce genre de manipulation reste très particulier
à chaque modèle considéré. Dans une approche automatique du problème, on
ne permettra qu'une analyse grossière. On regardera donc pour chaque
équation i, si celle-ci est donnée écrite de la forme x_i = second membre
fonction des autres variables (à laquelle on peut d'ailleurs se ramener
facilement si la variable est définie par son taux de croissance
(x_i / (t-1)) ou par son accroissement (x_i - x_{i-1}) ou encore par son logarithme).
D'après l'expérience acquise sur les modèles déjà étudiés, ceci semble
suffisant en première approche (voir le renvoi de la page 6).

La première partie de l'algorithmme (étapes 1 à 5) détermine le
prologue et l'épilogue du modèle. Elle se rapproche des études de Drud
[3] et Hellerman et Rarick [20]. La deuxième partie (étapes 6 à 9) qui
détermine les variables de bouclage est entièrement nouvelle.

L'algorithme :

### Etape 1 :

On détermine RC et CC qui sont respectivement les nombres d'éléments non nuls dans chaque ligne et chaque colonne de la matrice d'incidence. On pose $k=1$.

### Etape 2 :

Si aucun des RC(i) pour lesquels $k_i = 1$, n'est égal à 1, on va en 3. Sinon, on supprime la $i^{ème}$ ligne et la $j^{ème}$ colonne et on rajoute la $j^{ème}$ équation à l'ensemble d'équations qui constitue déjà le prologue comme la $k_i^{ème}$ élément de cet ensemble. On pose $k = k + 1$. On détermine à nouveau RC, CC et K (on rappelle que la $i^{ème}$ ligne et la $j^{ème}$ colonne ont été supprimées) et on retourne en 2.

### Etape 3 :

Si $k=n$, la matrice peut être mise sous la forme d'une matrice triangulaire inférieure en réordonnant les équations dans l'ordre où on a constitué le prologue qui est alors identique au système tout entier. Ceci signifie que le système est récursif et l'étude est terminée.

Sinon, on pose $k=n$ et on va en 4.

### Etape 4 :

Si aucun des CC(j), pour lesquels $k_j = 1$, n'est égal à un, on va en 5. Sinon, on supprime la $j^{ème}$ ligne et la $k_j^{ème}$ colonne, et on rajoute la $k_j^{ème}$ équation à l'ensemble des équations (numérotées de n à 1) qui constituent déjà l'épilogue comme $k_j^{ème}$ élément de cet ensemble. On pose $k = k-1$. On détermine à nouveau RC, CC et K et on retourne en 4.

### Etape 5 :

On réordonne les équations en plaçant au début les équations du prologue et à la fin les équations de l'épilogue dans l'ordre (indiqué par k) dans lequel on les a respectivement déterminées.

On détermine à nouveau K.

On ne considérera plus dans la suite que la matrice réduite correspondant au système constitué par les équations restantes.

Il s'agit maintenant de faire apparaître sur cette matrice les équations de bouclage. Le problème est donc de ramener la matrice qu'on vient de déterminer à la forme suivante (avec les notations de la figure 1) :
<figure><img src="image_5.png" /><figcaption>- **Figure 5** -</figcaption></figure>

où le nombre s d'équations de bouclage est le plus petit possible.

### Étape 6 :

La matrice d'incidence étant une matrice creuse (cf. Table 3), il est clair que le nombre s des variables de bouclage déterminées d'après la définition du paragraphe II.2 est en général faible par rapport à la dimension globale de la matrice. En réordonnant les équations de façon à placer à la fin du modèle les s équations définissant ces variables de bouclage, on peut ramener la matrice d'incidence à la forme décrite figure 5. Simultanément on redétermine K. Mais, a priori, cet ensemble de variables de bouclage qu'on vient de déterminer n'est pas minimal et un certain nombre de ces variables ne constitue que de "fausses" variables de bouclage. On teste donc chacune de ces variables (numérotées de 1 à s) de la façon suivante :

On dira que la variable $x_i$ est un prédecesseur de la variable $x_j$ si $i < j$ et s'il existe un certain nombre de variables $x_{h_1}, \dots, x_{h_p}$ avec
$$i < h_1 < h_2 < \dots < h_p < j$$

de telle sorte que les éléments $(h_1,i)$, $(h_2,h_1)$, \dots, $(j,h_p)$ de la matrice d'incidence ne soient pas nuls. (cf. figure 6).
<figure><img src="image_6.png" /><figcaption>- Figure 6 -</figcaption></figure>

D'autre part, on dira que la variable $x_i$ est un successeur de la variable $x_j$ si $i < j$ et si l'élément $(i,j)$ de la matrice d'incidence est non nul (cf. Figure 7).

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

On pose alors $k = 1$.

### Étape 7 :

Soit $x_t^i$ la variable de bouclage numérotée k.

Si $k_i = 0$, $x_t^i$ est une "vraie" variable de bouclage.

On supprime la ligne et la colonne correspondantes dans la matrice d'incidence. La suite de l'étude se fait sur la matrice réduite. On
détermine à nouveau K. On pose k = k+1 et on retourne en 7 sauf si toutes les variables de bouclage ont été testées, auquel cas l'étude est terminée. Sinon, on détermine l'ensemble des prédécesseurs et l'ensemble des successeurs de $x_t^i$. Si ces deux ensembles ont une intersection vide, on va en 8, sinon on va en 9.

### Etape 8 :

La variable de bouclage $x_t^i$ est une "fausse" variable de bouclage. On peut la faire disparaître en tant que variable de bouclage en réordonnant les équations de façon à avoir dans l'ordre l'ensemble des équations définissant l'ensemble des prédécesseurs de $x_t^i$ (ordre conservé), l'équation i et enfin l'ensemble des équations correspondant à l'ensemble des successeurs de $x_t^i$ (ordre conservé).

Dans l'exemple très simple suivant, où les variables sont numérotées de 1 à 4 :

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

l'ensemble des prédécesseurs de la variable 4 est ici l'ensemble $\{1,3\}$, l'ensemble des successeurs de la même variable est l'ensemble $\{2\}$. Ces deux ensembles sont distincts et en écrivant les équations dans l'ordre 1, 3, puis 4, puis 2, on obtient la nouvelle matrice :

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

Le résultat essentiel est que la variable $x_t^i$ qu'on a considéré, devient un prédecesseur de tous ses anciens successeurs ce qui supprime tous les 1 dans la matrice d'incidence au-dessus de l'élément de la
diagonale principale correspondant à cette variable. De plus, l'ordre étant conservé dans chacun des sous-ensembles qu'on a déterminé, on ne rajoute pas de 1 au-dessus de la diagonale principale.

On détermine à nouveau K.

On pose $s = s-1$, $k = k+1$ et on va en 7 sauf si toutes les variables de bouclage ont été testées, auquel cas l'étude est terminée.

### Étape 9 :

La variable de bouclage $x_t^i$ est une "vraie" variable de bouclage. Ceci correspond sur le graphe associé à l'existence d'une boucle sur cette variable. On supprime la ligne et la colonne correspondantes dans la matrice d'incidence. La suite de l'étude se fait sur la matrice réduite. On détermine à nouveau K. On pose $k = k+1$ et on va en 7, sauf si toutes les variables de bouclage ont été testées, auquel cas l'étude est terminée.

Dans l'étude du modèle de FAIR, par exemple, sur 83 variables, le nombre des variables de bouclage a été réduit de 30 à 6. Le temps de calcul global pour le programme de renumérotation écrit en FORTRAN H sur une machine IBM 370/168 a été de 1,30 seconde et laisse prévoir des temps de calcul tout à fait raisonnables pour des modèles de plus grande taille.

A titre d'exemple, on peut très aisément illustrer à la main les différentes étapes de l'algorithme proposé en l'appliquant au sous- modèle C.S [22] (conjectural-structural).

Ce petit modèle comporte 7 équations :

$$1. \quad S = 0.45 * A * (1+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$$

où les variables sont respectivement :
S : variation de salaire disponible

A : emploi salarié de l'industrie

L : taux de salaire horaire

P : prix à la consommation

H : productivité de la main d'oeuvre dans l'industrie

B : emploi dans l'industrie

D : variation de l'emploi dans l'industrie.

Le graphe associé à ce modèle est représenté sur la figure 8.
(pour les variables numérotées de 1 à 7) :

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

La matrice d'incidence du modèle donné est la suivante :

<table><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>RC</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>1</td><td>1</td><td></td><td></td><td></td><td></td><td>3</td></tr><tr><th>2</th><td></td><td>1</td><td></td><td></td><td></td><td></td><td>1</td><td>2</td></tr><tr><th>3</th><td></td><td></td><td>1</td><td>1</td><td>1</td><td></td><td></td><td>3</td></tr><tr><th>4</th><td></td><td></td><td>1</td><td>1</td><td>1</td><td></td><td></td><td>3</td></tr><tr><th>5</th><td></td><td></td><td></td><td></td><td>1</td><td>1</td><td></td><td>2</td></tr><tr><th>6</th><td></td><td></td><td></td><td></td><td></td><td>1</td><td>1</td><td>2</td></tr><tr><th>7</th><td></td><td></td><td>1</td><td></td><td></td><td></td><td>1</td><td>2</td></tr></tbody></table>

<table><caption>- Figure 9 -</caption><tr><th>CC</th><td>1</td><td>2</td><td>4</td><td>2</td><td>3</td><td>2</td><td>3</td></tr><tr><th>CC</th><td></td><td>1</td><td>3</td><td>2</td><td>3</td><td>2</td><td>3</td></tr><tr><th>CC</th><td></td><td></td><td>3</td><td>2</td><td>3</td><td>2</td><td>2</td></tr><tr><th>K</th><td></td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr></table>
On a reproduit figure 9 les différents calculs effectués sur RC et CC au cours des étapes 1 à 5. Tous les $k_i$ sont ici identiquement égaux à 1. La figure 10 représente le résultat obtenu à l'étape 5.

<figure><img src="image_11.png" /><figcaption>- figure 10 -</figcaption></figure>

On a pu ainsi déterminer le prologue qui est vide et l'épilogue qui est constitué par les équations 1 et 2.

Les étapes 6 à 9 sont décrites sur le système réduit dont la matrice d'incidence est la suivante :

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

D'après la définition du paragraphe II.2, les variables 4, 5, 6 et 7 sont des variables de bouclage. En effet, la valeur des $k_i$ n'intervient pas ici puisqu'ils sont tous égaux à 1, mais pour chacune des variables 4, 5, 6 et 7 la demi-colonne supérieure correspondante dans la matrice d'incidence contient des 1.

Par suite, l'ordre des équations et la matrice sont inchangés après l'étape 6.
Au cours des étapes 7 à 9, on teste successivement les variables 4, 5, 6 et 7.

Ainsi, pour la variable 4, l'ensemble des prédécesseurs est constitué par la variable 3 ; l'ensemble des successeurs est identique et la variable 4 est une "vraie" variable de bouclage. La matrice réduite devient :

<table><thead><tr><th></th><th>3</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><th>3</th><td>1</td><td>1</td><td></td><td></td></tr><tr><th>5</th><td></td><td>1</td><td>1</td><td></td></tr><tr><th>6</th><td></td><td></td><td>1</td><td>1</td></tr><tr><th>7</th><td>1</td><td></td><td></td><td>1</td></tr></tbody></table>

Pour la variable 5, l'ensemble des prédécesseurs est vide ; l'ensemble des successeurs est constitué par la variable 3 et la variable 5 est une "fausse" variable de bouclage. On peut la faire disparaître en tant que variable de bouclage en réordonnant les équations dans l'ordre : 5, puis 3, puis 6 et 7. La matrice d'incidence devient :

<table><thead><tr><th></th><th>5</th><th>3</th><th>6</th><th>7</th></tr></thead><tbody><tr><th>5</th><td>1</td><td></td><td>1</td><td></td></tr><tr><th>3</th><td></td><td>1</td><td></td><td></td></tr><tr><th>6</th><td></td><td></td><td>1</td><td>1</td></tr><tr><th>7</th><td></td><td>1</td><td></td><td>1</td></tr></tbody></table>

De même, la variable 6 est une "fausse" variable de bouclage. On réordonne les équations dans l'ordre 6, puis 5, puis 3 et 7, et la matrice d'incidence devient :

<table><thead><tr><th></th><th>6</th><th>5</th><th>3</th><th>7</th></tr></thead><tbody><tr><th>6</th><td>1</td><td></td><td></td><td>1</td></tr><tr><th>5</th><td></td><td>1</td><td></td><td></td></tr><tr><th>3</th><td></td><td></td><td>1</td><td></td></tr><tr><th>7</th><td></td><td></td><td></td><td>1</td></tr></tbody></table>
Pour la variable 7, l'ensemble des prédécesseurs est constitué par les variables 6, 5 et 3 ; l'ensemble des successeurs est constitué par la variable 6 et la variable 7 est une vraie variable de bouclage (les deux ensembles ont la variable 6 en commun).

Toutes les variables de bouclage ont été testées et l'étude est terminée. Il ne reste plus qu'à réordonner l'ensemble global des équations dans l'ordre dans lequel on a déterminé :

- le prologue : $\Phi$

- le coeur : équations 6, 5 et 3

- les équations de bouclage : équations 7 et 4

- l'épilogue : équations 2 et 1

La matrice d'incidence du système global devient :

<table><thead><tr><th></th><th>6</th><th>5</th><th>3</th><th>7</th><th>4</th><th>2</th><th>1</th></tr></thead><tbody><tr><th>Coeur</th><td>1</td><td></td><td></td><td>1</td><td></td><td></td><td></td></tr><tr><th></th><td>1</td><td>1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><th></th><td></td><td>1</td><td>1</td><td></td><td>1</td><td></td><td></td></tr><tr><th>Equation de bouclage</th><td>1</td><td></td><td></td><td>1</td><td></td><td></td><td></td></tr><tr><th></th><td>1</td><td>1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>Epilogue</th><td>1</td><td></td><td></td><td></td><td>1</td><td>1</td><td></td></tr><tr><th></th><td></td><td></td><td>1</td><td></td><td></td><td>1</td><td>1</td></tr></tbody></table>

matrice dont on rapprochera la structure de celle indiquée à la figure 1.
### III. RESOLUTION DU MODELE.

Il s'agit de résoudre, pour chaque période $t=1,...,T$, le système (4), qui est de dimension n.

Notre contribution à la résolution de ce problème est la mise au point d'une technique, proposée en 1976 (cf. [23]) et que nous décrivons au paragraphe III.1, qui permet de se ramener à la résolution d'un système de dimension s, s étant le nombre de variables de bouclage déterminé au chapitre précédent et l'évaluation des équations du système nécessitant la résolution d'un système récursif de dimension n-s. Cette technique est particulièrement efficace lorsque le rapport s/n est petit, et l'on a vu au chapitre II que c'est généralement le cas, en particulier pour le modèle STAR ($s/n = 0,02$).

Après avoir décrit en III.1 cette technique de réduction et avoir défini en III.2 quelques notations, nous décrirons plusieurs méthodes générales de résolution du système $f(x) = 0$: les méthodes de Newton, généralement considérées comme les plus efficaces, les méthodes de Gauss-Seidel et des moindres carrés qui, étant les plus simples à utiliser, sont les plus connues des économistes et la méthode des sécantes, qui a donné des résultats prometteurs en application au modèle STAR. Cette liste ne prétend pas être exhaustive. Nous avons, en particulier, laissé de coté l'algorithme de Brown [24], qui est une variante de la méthode de Newton basée sur l'élimination de Gauss et à convergence quadratique. L'efficacité de cette méthode sera étudiée par Rachdi [25].

Pour résoudre le système (4), on peut soit appliquer l'une de ces méthodes au système réduit de dimension s, et c'est ce que nous proposons, soit les appliquer directement au système de dimension n. Dans ce dernier cas, et si la méthode de Newton est retenue, il faut à chaque itération résoudre un système linéaire de dimension n. Dans les comparaisons numériques données à la fin de ce chapitre, nous avons utilisé un programme général (SSP-IBM) pour la résolution de ce système linéaire. Il est certain qu'il aurait été plus efficace d'utiliser un algorithme de résolution de système linéaire adapté aux matrices creuses, comme celui de Drud [3], mais nous n'en avons pas encore à notre disposition.
### III.1 Réduction de la taille du système.

Par définition des variables de bouclage et par construction des sous-systèmes (7), (8), (9) et (10), la situation à la période t est la suivante. Le vecteur $d_t$, qui est défini en (3), étant connu, on commence par résoudre le prologue (7) qui est récursif. Ensuite, pour une valeur donnée du vecteur $y_t = \{x_t^{m+1}, \dots, x_t^{m+s}\}$ des variables de bouclage, on sait que le coeur (8) du modèle est récursif et que l'on peut donc calculer $x_t^i$ pour $i = l+1, \dots, m$ par simple application successive des formules (8). Finalement, on connaît $x_t^i$ pour $i = 1, \dots, m+s$ et l'on peut alors évaluer, pour $i = m+1, \dots, m+s$, l'erreur sur les équations de bouclage (9), erreur égale à $k_i x_t^l - f_t^i(x_t^i, \dots, x_t^{m+s}, d_t)$.

Il est clair qu'une condition nécessaire et suffisante pour que $\{x_t^1, \dots, x_t^{m+s}\}$ soit solution de (8) - (9) est que cette erreur sur les équations de bouclage soit nulle. Puisque pour tout $y_t$ cette erreur peut-être évaluée récursivement par application successive des formules (8) et (9), nous considérerons que cette erreur est une fonction de $y_t$, notée $\varphi_t$. Le problème se ramène donc à la résolution du système de dimension $s$ :

$$ (11) \qquad \varphi_t^i(y_t) = 0, \quad i=1, \dots, s $$

On a vu que pour toute valeur numérique du vecteur $y_t$ on peut, à l'aide de (8) - (9), calculer la valeur numérique du vecteur $\varphi_t(y_t)$. Par contre, l'expression analytique des fonctions $\varphi_t^i$ est inaccessible, car pour l'obtenir, il faudrait effectuer l'élimination analytique des variables endogènes $x_t^i$, $i \le m$, entre les équations (8) et (9), ce qui est irréalisable dès que la taille du modèle dépasse la dizaine d'équations.

Nous utiliserons pour résoudre (11) une des méthodes itératives décrites ci-dessous. A chaque itération, la $k$-ième approximation $y_t^k$ étant connue, le calcul de $\varphi_t^k(y_t^k)$ se ramène à la résolution du système triangulaire (8) - (9), le prologue ayant été résolu avant le début des itérations. Une fois (11) résolu, on connaît $x_t^i$ pour $i=1, \dots, m+s$ et il ne reste plus, pour achever la résolution, qu'à résoudre l'épilogue (10) qui est, par construction, un système triangulaire.
Parmi les méthodes de résolution de (11) décrites ci-dessous,
deux, la méthode de Newton et celle des moindres carrés, font intervenir
les dérivées de $\varphi_t^i$, qui sont donc des dérivées de fonctions dont
l'expression analytique est inconnue. Nous allons décrire deux méthodes
pour le calcul de ces dérivées.

La première, la plus simple, consiste à approcher chaque élément
$\partial_t \varphi^i / \partial y_t^j$ de la matrice de Jacobi $\partial \varphi / \partial y$ du système (11) par une
différence finie :

$$
(12) \quad \frac{\partial \varphi_t^i}{\partial y_t^j} \approx \frac{1}{h} \left[ \varphi_t^i (y_t^1, \dots, y_t^{j-1}, y_t^{j+h}, y_t^{j+1}, \dots, y_t^s) - \varphi_t^i (y_t) \right]
$$

Pour chaque j, la résolution du système triangulaire (8) - (9) donne
une approximation de la $j$-ième colonne de la matrice $\partial\varphi / \partial y$. Par conséquent,
il faut résoudre $s+1$ fois le système (8) - (9) pour calculer la matrice
de Jacobi $\partial\varphi / \partial y$.

La deuxième méthode est une méthode exacte qui suppose connue
l'expression analytique des dérivées des équations du modèle. Nous allons
la détailler car l'expression (16) obtenue pour $\partial\varphi / \varphi \partial y$ sera à rapprocher
de l'expression du jacobien introduit dans la résolution du système
adjoint au chapitre IV.

La dérivation analytique des équations de (8) et (9) fait inter-
venir successivement :

- la matrice $A = ((\partial f_t^i / \partial x_t^j))$, $i,j = l+1,...,m$, dérivée des équations du coeur (8) par rapport aux variables du coeur ;

- la matrice $B = ((\partial f_t^i / \partial x_t^j))$, $i = l+1,...,m$ , $j = m+1,...,m+s$,
  dérivée des équations du coeur (8) par rapport aux variables de bouclage ;

- la matrice $C = ((\partial f_t^i / \partial x_t^j))$, $i = m+1,...,m+s$ , $j = l+1,...,m$,
  dérivée des équations de bouclage (9) par rapport aux variables du coeur ;

- la matrice $D = ((\partial f_t^i / \partial x_t^j))$, $i,j = m+1,...,m+s$, dérivée des équations de bouclage (9) par rapport aux variables de bouclage ;

- la matrice $\partial x / \partial y = ((\partial x_t^i / \partial y_t^j))$, $i = l+1,...,m$, $j = 1,...,s$, dérivée des variables du coeur par rapport aux variables de bouclage.
Avec ces notations, on déduit de (8) :

$$ (13) \qquad \frac{\partial x}{\partial y} = A \frac{\partial x}{\partial y} + B $$

Par ailleurs, avec $y_t = \{x_t^{m+1}, \dots, x_t^{m+s}\}$, on sait que les équations de bouclage s'écrivent :

$$ (14) \qquad \varphi_t^i = k_{m+i} y_t^i - f_t^{m+i} (x_t^1, \dots, x_t^m, y_t, e_t) , \quad i = 1, \dots, s $$

Avec les notations ci-dessus, on déduit immédiatement :

$$ (15) \qquad \frac{\partial \varphi}{\partial y} = L - C \frac{\partial x}{\partial y} - D $$

où L est la matrice diagonale de dimension s dont le ième élément diagonal est $k_{m+i}$. Eliminant $\partial x / \partial y$ entre (13) et (15), on obtient l'expression analytique du jacobien $\partial \varphi / \partial y$ du système (11) :

$$ (16) \qquad \frac{\partial \varphi}{\partial y} = L - C(I-A)^{-1}B - D $$

Pour calculer $\partial \varphi / \partial y$, on commence par résoudre le système linéaire (13), ce qui est simple, puisque l'on sait d'après (8) que la matrice A est triangulaire inférieure. Connaissant $\partial x / \partial y$, le jacobien est donné par application de la formule (15). La seule difficulté est donc le calcul _analytique_ que représente la dérivation des équations du modèle. Pour des modèles importants, il est irréaliste de faire ce calcul à la main, c'est pourquoi nous utilisons actuellement la formule approchée (12). Nous projetons cependant, dans le cadre du projet MODULECO [1], de construire un préprocesseur qui, à partir des expressions FORTRAN de calcul des expressions $k_i x_t^i - f_t^i (x_t, e_t)$ données par l'utilisateur, génèrerait les expressions FORTRAN permettant le calcul de ces dérivées. Ce préprocesseur sera construit en reprenant les bonnes idées des dérivateurs existants (FORMAC et celui du logiciel XING) et en tenant compte du fait que le nombre de fonctions apparaissant dans un modèle macroéconomique est très faible (quatre opérations arithmétiques, logarithmes et exponentielles).
### III.2. Notations.

Les méthodes de résolution d'un système algébrique non linéaire qui seront décrites aux paragraphes suivants peuvent être aussi bien appliquées au système global (4) qu'au système réduit (11). Nous adopterons par conséquent une terminologie générale inspirée de [26].

Soit $f$ une application de $\mathbb{R}^n$ dans $\mathbb{R}^n$. Il s'agit de résoudre le problème :

$$ (17) \qquad f(x) = 0 $$

Nous noterons $x^*$ toute solution du problème(17), i.e. tout $x$ tel que $f(x) = 0$. Nous poserons :

$$ (18) \qquad f'(x) = \left( \frac{\partial f_i(x)}{\partial x_j} \right) ; \quad f''(x) = \left( \frac{\partial^2 f_i(x)}{\partial x_j \partial x_k} \right) $$

Nous utiliserons les normes suivantes :

$$ (19) \qquad \begin{gathered} \|x\| = \max_i |x_i| \quad ; \quad \|f'(x)\| = \max_i \sum_{j=1}^n \left| \frac{\partial f_i(x)}{\partial x_j} \right| \\ \|f''(x)\| = \max_{i,j} \sum_{k=1}^n \left| \frac{\partial^2 f_i(x)}{\partial x_j \partial x_k} \right| \end{gathered} $$

Pour tout $x \in \mathbb{R}^n$ et tout scalaire positif $a$, nous appelerons $V(x,a)$ le pavé fermé de centre $x$ et de rayon $a$ :

$$ (20) \qquad V(x,a) = \{y \in \mathbb{R}^n / \|x-y\| \le a\} = \{y \in \mathbb{R}^n / |x_i - y_i| \le a, i=1, \dots, n\} $$

Enfin, les méthodes décrites ci-dessous étant toutes des méthodes itératives, nous noterons $x^0$, $x^k$ l'approximation initiale, $x^{k+1}$ l'approximation connue au début de l'itération $k$ et $e^k$ l'erreur :

$$ (21) \qquad e^k = x^k - x^* $$

### III.3. La méthode de Newton.

A la $k^{\text{ième}}$ itération, connaissant $x^k$, on calcule $x^{k+1}$ en résolvant le système linéaire suivant :

$$ (22) \qquad f'(x^k)(x^{k+1} - x^k) = -f(x^k) $$

Nous admettrons le résultat théorique suivante (cf. [26]) :
**Théorème 2 :**

Soit $x^0$ l'approximation initiale. On suppose que $f'(x^0)$ est
inversible et l'on note :

$$
(23) \quad a = \|[f'(x^0)]^{-1}\| ; \quad b = \|[f'(x^0)]^{-1}f(x^0)\| ; \quad W = V(x^0, 2b)
$$

On suppose que $f$ est de classe $C^2$ dans le domaine $W$ ; on définit
les scalaires $c$ et $\mu$ par :

$$
(24) \quad c = \sup_{x \in W} \|f''(x)\| ; \quad \mu = 2n a b c
$$

Alors, $\mu \le 1$ est une condition <u>suffisante</u> pour que la suite $\{x^k\}$
obtenue par l'algorithme (22) soit contenue dans $W$ et converge vers
un $x^*$ qui est solution de $f(x) = 0$. On a une majoration de l'erreur :

$$
(25) \quad \| e_k \| \le \left(\frac{1}{2}\right)^{k-1} \mu^{2k-1} b , \forall k \ge 0
$$

Si $\mu < 1$, on dit que la convergence est "superrapide".

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

Ce résultat est, à notre connaissance, le plus général, car c'est
le seul qui ne fasse aucune hypothèse sur l'inversibilité de f'(x)
pour x = x* , hypothèse gênante puisque x* est inconnu... Le théorème 2
a cependant des défauts :

1) les scalaires a, et surtout c, sont difficiles à calculer ;

2) la condition de convergence $\mu = 2n a b c \le 1$ est très sévère.
Par exemple, avec $f(x) = x^3 - 1$ et $x^0 = 4$, on a $\mu = 1,8$ alors
que l'algorithme (22) converge.

3) le théorème 2 n'assure pas la convergence **quadratique** de la méthode. Par exemple, avec $f(x) = x^2$ et $x^0 = 1$, on a $\mu = 1$, le théorème s'applique, mais la convergence n'est que du premier ordre.

Pour garantir la convergence quadratique de la méthode de Newton,
il faut faire des hypothèses sur l'inversibilité du jacobien au voisi-
nage de la solution :

Théorème 3 :

Supposons a priori que le problème $f(x) = 0$ ait une solution $x^*$.
Soit $x^0$ l'approximation initiale et soient $\alpha, \beta$ et $\gamma$ les scalaires
définis par :
$$
(26) \quad \| [f'(x)]^{-1} \| \le \alpha \quad \| f''(x) \| \le \frac{2\beta}{n}, \forall x \in V(x^*, \|e^0\|)
$$

$$
(27) \qquad \gamma = \alpha \beta \|e^0\|^2
$$

Alors $\gamma < 1$ est une condition suffisante pour que la suite $\{x^k\}$ converge vers $x^*$ et, si elle est vérifiée, la convergence est quadratique avec :

$$
(28) \quad \| e^{k+1} \| \le \alpha \beta \| e^k \|^2 , \forall k \ge 0 .
$$

Le théorème s'applique si l'approximation initiale $x^0$ appartient
à un certain voisinage de la solution $x^*$, c'est pour cela que l'on
parle de **convergence localement quadratique**.

Pour utiliser la méthode, on peut se servir directement de l'algo-
rithme (22), en résolvant un système linéaire de dimension n à chaque
itération, ou bien remplacer (22) par l'algorithme suivant, qui lui est
équivalent :

$$
(29) \quad x^{k+1} = x^k - [f'(x^k)]^{-1} f(x^k)
$$

où la matrice $[f'(x^k)]^{-1}$ est calculée à l'aide de la connaissance de
$[f'(x^{k-1})]^{-1}$ par une méthode proposée par Powell [27].

### III.4. La méthode de Newton modifiée.

L'application de l'algorithme (22) exige le calcul du jacobien $f'(x^k)$ à chaque itération. Ce calcul pouvant être long, il est naturel lorsque l'on sait que $x^0$ est proche de la solution $^{(1)}$, de penser à remplacer, dans (22), $f'(x^k)$ par $f'(x^0)$ : on ne calculera le jacobien que pour la première itération. L'algorithme est alors le suivant :

$$
(30) \quad f'(x^0)(x^{k+1} - x^k) = -f(x^k)
$$

<footnote>(1) ce qui est notamment le cas lorsque l'on résout un modèle macro-
économique pour un période : l'approximation initiale est la
solution de la période précédente.</footnote>
On peut, naturellement, calculer une fois pour toutes l'inverse
[f'(x°)]⁻¹ du jacobien mais, en général, le nombre d'itérations à
effectuer est trop petit pour que cette opération soit rentable. On
démontre :

Théorème 2' : (cf. [26]) :

Avec les notations du théorème 2, μ < 1 est une condition suffi-
sante pour que l'algorithme (30) converge vers la solution de f(x) = 0
avec l'estimation suivante de l'erreur :

$$
(31) \qquad e^k \le 2b \mu^k
$$

Théorème 3' :

Avec les notations et les hypothèses du théorème 3, l'algorithme
(30) converge vers la solution de f(x) = 0 avec l'estimation suivante
de l'erreur :

$$
(32) \quad e^{k+1} \le \gamma \|e^k\|
$$

On notera que α et β ne dépendant que de la régularité du problème,
γ étant égal à αβ||e°||, plus x° est proche de x*, plus γ est petit et
plus la convergence (32) est rapide.

### III.5. La méthode de Gauss-Seidel.

L'algorithme :

$$
(33) \quad x^{k+1} = x^k + a f(x^k)
$$

de résolution de $f(x) = 0$ en dimension 1, où a est un paramètre, est
bien connu. Avec les notations de III.2, on définit q par :

$$
(34) \quad |1 + af'(x)| \le q, \quad \forall x \in V(x^*, |e^0|)
$$

et l'on a le résultat classique :

$$
(35) \quad |e^{k+1}| \le q |e^k|, \quad \forall k \ge 0
$$

qui montre que q < 1 est une condition suffisante de convergence de (33)
et que celle-ci, lorsqu'elle a lieu, est du premier ordre.
La méthode de Gauss-Seidel est une généralisation de l'algorithme (33) au cas multidimensionnel. Cette généralisation donne naissance à deux algorithmes distincts :

$$ (36) \qquad x_i^{k+1} = x_i^k + a_i f_i(x_i^k), \quad i = 1, \dots, n $$

$$ (37) \qquad x_i^{k+1} = x_i^k + a_i f_i(x_1^{k+1}, \dots, x_{i-1}^{k+1}, x_i^k, \dots, x_n^k), \quad i = 1, \dots, n $$
où les $a_i$ sont des paramètres.

L'algorithme (36) est la généralisation naturelle de (33) mais il n'est jamais utilisé : il est moins efficace que (37) et, contrairement à l'apparence, il est aussi plus difficile à programmer.

Un troisième algorithme, "type Gauss-Seidel", est le suivant :

$$ (38) \qquad f_i(x_1^{k+1}, \dots, x_i^{k+1}, x_{i+1}^k, \dots, x_n^k) = 0 \quad i = 1, \dots, n $$

A chaque itération et pour chaque $i$, connaissant $x_1^{k+1}, \dots, x_{i-1}^{k+1}$ d'une part et $x_{i+1}^k, \dots, x_n^k$ d'autre part, on résout l'équation n° $i$ par rapport à la variable n° $i$. Une itération revient donc à la résolution de n problèmes $f(x) = 0$ en dimension $l$. Cette méthode est surtout intéressante lorsque la ième équation est analytiquement résoluble en $x_i$. En général, ceci est vrai pour certaines équations et pas pour d'autres et "la" méthode de Gauss-Seidel généralement utilisée est une combinaison des algorithmes (37) et (38). Soit $I$ l'ensemble des indices des équations résolubles par rapport à la variable portant le même numéro. Pour tout $i \in I$, soit $h_i$ la fonction telle que $f_i(x) = 0$ soit équivalent à $x_i = h_i(x)$. Le meilleur algorithme de Gauss-Seidel à notre connaissance, est alors le suivant :

$$ (39) \qquad \left\{ \begin{array}{ll} x_i^{k+1} = h_i(x_1^{k+1}, \dots, x_{i-1}^{k+1}, x_{i+1}^k, \dots, x_n^k), & \text{si } i \in I \\ \\ x_i^{k+1} = x_i^k + a_i f_i(x_1^{k+1}, \dots, x_{i-1}^{k+1}, x_i^k, \dots, x_n^k), & \text{si } i \notin I \end{array} \right\} \quad i=1, \dots, n $$

L'étude de la convergence de ces algorithmes fait intervenir des conditions sur les valeurs propres de la matrice $f'(x)$ dans le voisinage de la solution $x^*$. En pratique donc, il est totalement illusoire de chercher à s'assurer de la convergence a priori de la méthode de Gauss-Seidel et l'ajustement, c'est à dire le choix des "poids" $a_i$, qui
joue un rôle fondamental pour la convergence, est fait "à tâtons", ce qui
rend la méthode de Gauss-Seidel, si simple en apparence, d'usage délicat.
Il est clair, par ailleurs, que la convergence de Gauss-Seidel ne peut
être que du premier ordre.

### III.6. La méthode des sécantes.

Il s'agit d'une méthode dont la convergence est plus rapide que
celle de Gauss-Seidel et qui, contrairement à Newton, ne fait pas
intervenir les dérivées de f.

Le principe de la méthode est facile à comprendre en dimension 1,
l'algorithme s'écrivant alors :

$$
(40) \qquad x^{k+1} = x^k - \frac{x^k - x^{k-1}}{f(x^k) - f(x^{k-1})} f(x^k)
$$

L'interprétation géométrique de l'algorithme (40) est simple :
l'approximation est l'abcisse de l'intersection avec Ox de la droite
joignant les deux points de la courbe y = f(x) d'abcisse x^k et x^{k-1}.

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

On constatera que lorsque $x^k$ est proche de sa limite $x^*$, les points $x^k$ et $x^{k-1}$ sont voisins et le quotient $(x^k - x^{k-1}) / [f(x^k) - f(x^{k-1})]$ est proche de $[f'(x^k)]^{-1}$ : l'algorithme (40) est alors proche de celui de Newton. On démontre aisément le :

Théorème 4 :

Supposons que la suite {$x^k$} soit contenue dans un intervalle [a,b]
et que la fonction f vérifie :

$$
(41) \quad |f'(x)| \ge \alpha, \quad |f''(x)| \le \beta, \quad \forall x \in [a,b]
$$
alors l'algorithme (40) vérifie :

$$
(42) \quad |e^{k+1}| \le \frac{\beta}{2\alpha} |e^k| |e^{k-1}|, \forall k \ge 0
$$

Ce résultat montre qu'au voisinage de $x^*$, la convergence est
quasi-quadratique.

Passons maintenant au cas général de la résolution de $f(x) = 0$ dans $\mathbb{R}^n$. La méthode est basée sur une généralisation de l'interprétation géométrique donnée ci-dessus, à savoir : connaissant $n+1$ approximations $x^{k-n}$, $x^{k-n+1}$, ..., $x^k$ de $x^*$, on remplace $f(x)$ par l'application linéaire (unique) $L(x) = Ax + b$ qui vérifie $L(x^{k-j}) = f(x^{k-j})$, $j=0, \dots, n$ et l'on choisit pour $x^{k+1}$ la solution du système linéaire $L(x^{k+1}) = 0$.

On montre que l'algorithme est alors le suivant :

$$
(43) \quad x^{k+1} = x^k - H_k G_k^{-1} f(x^k), \quad \forall k \ge n
$$

avec les notations :

$$
(44) \quad \begin{gathered} \{H_k\}_{ij} = x_i^{k-j+1} - x_i^{k-j} \quad , \quad \{G_k\}_{ij} = f_i(x^{k-j+1}) - f_i(x^{k-j}) \\ i,j = 1, \dots, n \end{gathered}
$$

Pour l'énoncé du résultat de convergence, qui est assez complexe,
ainsi que pour la démonstration, nous renvoyons à Rachdi [25] qui, en
plus d'un certain nombre d'hypothèses de régularité sur f, suppose que
les matrices initiales G_n et H_n sont inversibles.

L'algorithme (43) pose un problème évident, c'est celui de
l'inversion de la matrice G_k. Nous proposons l'algorithme suivant pour
calculer G_{k+1}^{-1} connaissant G_k^{-1} (cf. [25] pour la démonstration). On
calcule d'abord le vecteur z = G_k^{-1} [f(x^{k+1}) - f(x^k)] et l'on en déduit
G_{k+1}^{-1} par :

$$
(45) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\{G_{k+1}^{-1}\}_{1,j} & = \frac{1}{z_n} \{G_k^{-1}\}_{n,j} \quad , \quad j=1,\dots,n \\[1em]
\{G_{k+1}^{-1}\}_{i,j} & = \{G_k^{-1}\}_{i-1,j} - \frac{z_{i-1}}{z_n} \{G_k^{-1}\}_{n,j} \quad , \quad i=2,\dots,n \quad , \quad j=1,\dots,n
\end{array}
\right.
$$
On peut montrer que $det(G_{k+1}) = z_n det(G_k)$. D'après (43), une condition nécessaire de convergence est l'inversibilité de toutes les matrices $G_k$, qui implique donc que $z_n$ n'est jamais nul et que la formule (45) est toujours applicable.

On notera que la formule (45) est moins coûteuse que celle de Ortega-Rheinboldt [28].

### III.7. Méthode des moindres carrés.

Soit $F(x)$ l'application de $\mathbb{R}^n$ dans $\mathbb{R}$ définie par :

$$ (46) \qquad F(x) = \sum_{i=1}^{n} f_{i}^{2}(x) $$

Il est clair que $F(x)$ est nul si et seulement si $f(x) = 0$, c'est à dire si $x$ est solution de notre problème. On notera que le gradient $F'(x)$ de $F$ se déduit du jacobien $f'(x)$ de $f$ par :

$$ (47) \qquad F'(x) = 2f'(x) f(x) $$

On en déduit le résultat suivant, dont la démonstration est immédiate :

### Théorème 5 :

On suppose que $f$ est différentiable de $\mathbb{R}^n$ dans $\mathbb{R}^n$. On a alors :

1) le problème de la minimisation sans contraintes de $F$ a toujours au moins une solution.

2) tout $x$ vérifiant $F(x) = 0$ est solution du système $f(x) = 0$.

3) si $\hat{x}$ est un minimum local de $F$ tel que $f'(\hat{x})$ est inversible, alors $\hat{x}$ est solution de $f(x) = 0$.

On voit donc que l'on retombe sur la même difficulté que pour la méthode de Newton, à savoir l'inversibilité du jacobien du système. Si à l'itéré $x^k$, le déterminant de $f'(x^k)$ est quasi-nul, alors Newton diverge et la méthode des moindres carrés "se bloque".

On utilise, pour minimiser $F$, la méthode du gradient conjugué ou celle de Davidon-Fletcher-Powell (cf. chapitre V). On notera que toutes deux exigent le calcul du gradient $F'(x)$, donc du jacobien $f'(x)$.
### III.8. Comparaisons numériques.

Les tests présentés ici ont tous été effectués sur le modèle STAR [7]. Des tests antérieurs [23] avaient été effectués sur le petit modèle PIMPON [6] mais ceux sur STAR sont plus significatifs. Le modèle STAR comporte 139 équations, mais les 48 premières équations ne contiennent pas le contrôle et forment un système récursif résolu une fois pour toutes au moment de la lecture des données. Les tests que nous présentons portent donc sur un modèle de 91 équations, que nous avons systématiquement résolu pour $T = 5$ périodes.

Pour comparer valablement les algorithmes de résolution, nous avons utilisé la procédure suivante. Nous résolvons $K$ fois le modèle pour $K$ valeurs différentes du vecteur contrôle $u$ ; le vecteur $u$ est à chaque fois déterminé de la manière suivante : chacune de ses composantes est l'observation d'une variable aléatoire uniformément répartie sur l'intervalle $[0,q]$, où $q$ est un paramètre. On choisit systématiquement comme approximation initiale la solution du modèle qui correspond à $u = 0$. Par conséquent, avec une petite valeur de $q$, l'approximation initiale est proche de la solution et avec une grande valeur de $q$, l'approximation initiale est mauvaise. En faisant varier le paramètre $q$, nous obtenons ainsi un temps moyen de résolution du modèle fonction de la qualité de l'approximation initiale.

Pour éviter les problèmes dûs aux facteurs d'échelle, nous avons effectué un changement systématique d'unités de façon à ce que toutes les variables endogènes du modèle STAR soient voisines de l'unité.

Pour les méthodes de résolution globales qui s'appliquent directement à la résolution du système (4), nous avons choisi pour chaque période le test d'arrêt suivant :

$$ (48) \qquad \sqrt{\frac{1}{n} \sum_{i=1}^{n} [k_i x_{it}^i - f_t^i(x_t, d_t)]^2} \le \varepsilon $$

Pour les méthodes utilisant la technique de réduction décrite au paragraphe III.1, il s'agit de résoudre le système (11) et le test d'arrêt pour chaque période est le suivant :
$$ (49) \qquad \sqrt{\frac{1}{s} \sum_{i=1}^{s} [\varphi_t^i(y_t)]^2} \le \varepsilon $$

Grâce à notre changement d'unités, nous sommes assurés que les tests (48) et (49) sont équivalents.

Pour les méthodes qui nécessitent le calcul du jacobien du système, ce dernier a été calculé par différences finies et, en particulier, les dérivées partielles du système (11) ont été approchées par la formule (12).

Enfin, précisons que toutes les méthodes ont été programmées en FORTRAN H et en double précision sur l'ordinateur IBM 370-168 du C.I.R.C.E.

La méthode de Newton "globale", c.a.d. appliquée directement au système (4), et programmée soit à l'aide de l'algorithme (22), soit en utilisant le programme NSO1AD de la bibliothèque HARWELL [29], qui est basée sur l'algorithme (29), a donné des résultats catastrophiques : dans les deux cas, un temps de résolution de l'ordre de 18 secondes à comparer avec les résultats de la table 5. Il est clair que ce résultat peut être amélioré si l'on tient explicitement compte que la matrice $f'(x^k)$ est creuse pour la calculer et pour résoudre le système linéaire (22).

L'initialisation de la méthode des sécantes exige le calcul de $G_n$ (n résolutions du système récursif) puis l'inversion de la matrice $G_n$. En application au système global (4), avec $n=91$, cette initialisation est très coûteuse et bien que l'application de la formule (45) rende la méthode assez rapide, initialisation non comprise, la méthode des sécantes "globale" converge en un temps moyen de 20 secondes, ce qui est également catastrophique.

Nous avons utilisé le programme écrit par les auteurs du modèle STAR, basé sur l'algorithme (39) qui est, à notre avis, la meilleure version de l'algorithme de Gauss-Seidel. Nous avons conservé les valeurs des paramètres $a_i^k$ de (39) qui ont été optimisées par les auteurs du modèle. Avec $q = 10^{-3}$, c.a.d. avec une bonne approximation initiale, on obtient pour différentes valeurs de la précision $\epsilon$ les temps de calcul suivants :
<table><thead><tr><td>ε</td><td>temps de calcul</td></tr></thead><tbody><tr><td>10<sup>-2</sup></td><td>0,097 sec.</td></tr><tr><td>10<sup>-4</sup></td><td>0,344 sec.</td></tr><tr><td>10<sup>-7</sup></td><td>0,762 sec.</td></tr></tbody></table>

### - Table 4 - Temps de résolution de Gauss-Seidel

Nous avons retenu cinq méthodes basées sur la technique des variables de bouclage, c.a.d. appliquées au système (11).

* Méthode M1 : méthode de Newton classique basée sur l'algorithme (22).

* Méthode M2 : programme NSO1AD de la bibliothèque HARWELL [27]: variante de Newton basée sur l'algorithme (29), le calcul de $[f'(x^k)]^{-1}$ étant effectué à l'aide de la connaissance de $[f'(x^{k-1})]^{-1}$ par l'algorithme de Broyden (cf. [29]).

* Méthode M3 : méthode de Newton modifiée basée sur l'algorithme (30).

* Méthode M4 : méthode des sécantes basée sur l'algorithme (43) avec utilisation de la formule (45) pour calculer $G_{k+1}^{-1}$ à partir de $G_k^{-1}$.

* Méthode M5 : programme VA05AD de la bibliothèque HARWELL : méthode des moindres carrés, c.a.d. minimisation du critère (46) par la méthode de Fletcher-Powell.

La table 5 donne les temps de résolution moyens en fonction de la précision $ε$ et de la qualité moyenne $q$ de l'approximation initiale :
<table><caption>- Table 5 - Comparaison de temps de résolution.</caption><thead><tr><th>ε</th><th>q</th><th>M1</th><th>M2</th><th>M3</th><th>M4</th><th>M5</th></tr></thead><tbody><tr><td rowspan="3">10<sup>-2</sup></td><td>10<sup>-3</sup></td><td>0,042 sec.</td><td>0,011 sec.</td><td>0,037 sec.</td><td>0,011 sec.</td><td>0,044 sec.</td></tr><tr><td>1</td><td>0,067 sec.</td><td>0,077 sec.</td><td>0,043 sec.</td><td>0,047 sec.</td><td>0,081 sec.</td></tr><tr><td>100</td><td>0,101 sec.</td><td>0,087 sec.</td><td>0,050 sec.</td><td>0,057 sec.</td><td>0,090 sec.</td></tr><tr><td rowspan="3">10<sup>-4</sup></td><td>10<sup>-3</sup></td><td>0,067 sec.</td><td>0,062 sec.</td><td>0,044 sec.</td><td>0,042 sec.</td><td>0,066 sec.</td></tr><tr><td>1</td><td>0,067 sec.</td><td>0,082 sec.</td><td>0,044 sec.</td><td>0,046 sec.</td><td>0,086 sec.</td></tr><tr><td>100</td><td>0,099 sec.</td><td>0,089 sec.</td><td>0,059 sec.</td><td>0,059 sec.</td><td>0,100 sec.</td></tr><tr><td rowspan="3">10<sup>-7</sup></td><td>10<sup>-3</sup></td><td>0,068 sec.</td><td>0,061 sec.</td><td>0,044 sec.</td><td>0,044 sec.</td><td>0,066 sec.</td></tr><tr><td>1</td><td>0,100 sec.</td><td>0,081 sec.</td><td>0,050 sec.</td><td>0,055 sec.</td><td>0,090 sec.</td></tr><tr><td>100</td><td>0,110 sec.</td><td>0,098 sec.</td><td>0,065 sec.</td><td>0,068 sec.</td><td>0,104 sec.</td></tr></tbody></table>
Deux méthodes se dégagent de notre test : M3 (Newton modifiée) et M4 (méthode des sécantes). La méthode de Newton modifiée, qui diverge parfois avec $q=1000$ (mauvaise approximation initiale), semble cependant moins sûre. Dans l'état actuel de nos connaissances, nous retiendrons comme meilleure méthode pour résoudre le modèle STAR la méthode des sécantes appliquée au système (11).

La table 4 ayant été obtenue pour $q=10^{-3}$, on constate que la méthode de Gauss-Seidel est battue par toutes les méthodes testées à la table 5 et ceci d'autant plus que la précision demandée est plus grande ($ε$ plus petit). Rappelons à ce sujet que, au moins si la résolution du modèle n'est qu'une partie d'un programme d'optimisation utilisant des méthodes sophistiquées (Fletcher-Powell, gradient conjugués, etc...), il est indispensable de résoudre le modèle avec grande précision pour que l'algorithme d'optimisation fonctionne normalement.

Rappelons enfin en conclusion que les méthodes testées à la table 5, lorsqu'elles sont appliquées directement au système global (4), donnent des temps de calcul supérieurs à 15 secondes, ce qui démontre de façon indiscutable l'intérêt de la technique des variables de bouclage.
# IV. OPTIMISATION DU MODELE.

## IV.1. Notations, définitions et validité de l'étude.

Les entiers n, p, r et q étant toujours ceux définis au début du paragraphe II.1, on pose $M = (p+1)n + (q+1)r$. Le critère, ou fonction objectif du modèle, ou fonction coût, est donné comme une application de $R^M$ dans $R$ :

$$ (50) \qquad j(x,u) = \sum_{t=1}^{T} j_t(x_t, x_{t-1}, \dots, x_{t-p}, u_t, u_{t-1}, \dots, u_{t-q}) $$

Le problème étudié dans ce chapitre est la minimisation du critère (50) compte tenu du fait que x et u sont liés par le système (1).

La formulation standard du problème est la suivante. On suppose que toutes les fonctions définissant le modèle sont continues, que pour tout $u \in R^R$ le modèle a une solution et une seule $x(u)$ et l'on minimise le critère $J(u) = j[x(u),u]$.

Ce cadre théorique est ici inapplicable. En effet, les fonctions décrivant le modèle contiennent souvent des logarithmes et toujours des quotients de variables : elles ne sont ni continues ni même définies pour toutes les valeurs possibles de leurs arguments. Il est de plus impossible d'affirmer a priori que la solution du modèle est unique.

Remarquons par ailleurs qu'un modèle macroéconomique n'est pas une loi absolue comme en physique, mais une estimation du comportement réel de l'économie dans un voisinage d'un état de référence de cette économie obtenu par extrapolation du passé. Il est bien connu qu'après un choc brutal (mai 68, crise du pétrole etc...), les modèles doivent être réestimés. Par contre, dans une période "stable", on peut estimer que le comportement de l'économie est dynamiquement continu : la solution $x_t$ à la période t est "proche" de la solution $x_{t-1}$ de la période précédente.

Nous nous limiterons donc volontairement au cadre théorique suivant. On suppose connue une politique de contrôle de référence $\bar{u}$. On limite la recherche du contrôle optimal à un voisinage $U$ de $\bar{u}$ que l'on suppose fermé et borné dans $R^R$. Pour tout $u \in U$, on suppose que le modèle (1) a au moins une solution. Pour tout $u \in U$, on appelle $x(u)$ la trajectoire définie de la manière suivante : $x_1(u)$ est la solution du système (1)
pour la période 1 la plus proche dans $\mathbb{R}^n$ de l'état initial connu $x_0(^*)$,
$x_2(u)$ est la solution pour $t=2$ la plus proche de $x_1(u)$, et ainsi de
suite. On définit ainsi une application de $U$ dans $X = x(U)$ et l'on pose :

$$
(51) \quad J(u) = j[x(u),u]
$$

et le problème fondamental est de minimiser J sur U :

Problème P : trouver $u^* \in U$ tel que :

$$
(52) \quad J(u^*) \le J(u), \quad \forall u \in U
$$

Les ensembles U et X étant des voisinages de la politique de
référence ū et de la trajectoire de référence x(ū), on peut supposer
que les fonctions f_t définissant le modèle (1) et les fonctions j_t
définissant le critère (50) sont continues pour tout couple (u,x) de
U.X . On démontre alors aisément (U est fermé, borné, non vide) que
le problème (52) admet au moins une solution.

Il faut noter que même si le critère (50) est quadratique, la
solution de (52) n'est pas forcément unique, vu la non-linéarité du
modèle, mais un algorithme itératif d'optimisation partant de la poli-
tique de référence ū convergera vers l'optimum local le plus proche de
ū, ce qui est suffisant en pratique.

Il existe deux façons d'aborder le problème (52). La première
méthode, que nous appellerons **méthode implicite**, consiste à minimiser
**directement** le critère $J(u)$ qui est défini **implicitly** par la donnée
du modèle. Lorsque la minimisation de $J$ s'effectue par une méthode de
gradient, cette méthode s'appelle la **méthode du gradient réduit** [30].
A chaque itération d'une méthode implicite, le modèle est résolu
exactement. La deuxième méthode est la **méthode de pénalisation** [31], qui
consiste à considérer les équations du modèle comme des contraintes et
à minimiser sous ces contraintes le critère $j(x,u)$ par rapport à $x$ et $u$,
le critère $j$ étant connu, lui, **explicitement**.

La méthode implicite la plus simple est la méthode séquentielle [32]

<footnote>(*) dans le cas bien extraordinaire où deux solutions du modèle seraient
à égale distance de $x_0$, on choisit au hasard.</footnote>
exposée au paragraphe IV.2. L'application de la méthode du gradient réduit nécessite l'utilisation d'un algorithme standard de minimisation sans contraintes d'une fonction de plusieurs variables et d'une méthode de calcul du gradient J'(u) de la fonction implicite J(u). Nous donnerons ci-dessous deux méthodes pour calculer ce gradient.

Nous allons donner quelques définitions et notations utiles dans les paragraphes ultérieurs.

Pour tout $u \in U$, le gradient $J'(u)$ de la fonction (51) est le vecteur de $R^R$ qui vérifie :

$$ (53) \qquad (J'(u),v) = \lim_{h \to 0} \frac{1}{h} [J(u+hv) - J(u)] , \forall v \in R^R, \\ \text{où } ( , ) \text{ est le produit scalaire de } R^R. $$

Soient $h_i$, $i=1,...,R$, des scalaires quelconques. Nous appellerons $\delta_i$ le vecteur de $R^R$, de composantes $\delta_{ij}$, avec :

$$ (54) \qquad \left\{ \begin{array}{ll} \delta_{ii} = h_i & i,j = 1,...,R \\ \delta_{ij} = 0, & \text{si } i \neq j \end{array} \right. $$

Pour tout vecteur-fonction $f$ d'un vecteur $x$ mais ne dépendant pas nécessairement uniquement de $x$, on appelle $\partial f / \partial x$ la matrice des $\partial f_i / \partial x_j$. Les dimensions de $f$ et $x$ sont quelconques, éventuellement égales à 1, donc $\partial f / \partial x$ désigne suivant les cas une matrice, un vecteur ou un scalaire.

On appelle $B$ l'ensemble des variables de bouclage du modèle, définies au chapitre II. Pour tout vecteur $z$ de $R^n$, on pose :

$$ (55) \qquad z^a = \{z_i\}_{i \in B} , \quad z^b = \{z_i\}_{i \in B} $$

Par définition des variables de bouclage, le modèle (1) s'écrit pour la période $t$ :

$$ (56) \qquad \left\{ \begin{aligned} x_t^a &= f_t^a(x_t, \dots, x_{t-p}, u_t, \dots, u_{t-q}) \\ Lx_t^b &= f_t^b(x_t, \dots, x_{t-p}, u_t, \dots, u_{t-q}) \end{aligned} \right. $$

où la matrice $\partial f_t^a / \partial x_t^a$ est **triangulaire inférieure** et où $L$ est la matrice diagonale de dimension $s$ d'éléments diagonaux $k_{i+m}^s$.
IV.2. La méthode séquentielle.

A chaque étape n°i de l'itération, on minimise la fonction
scalaire d'une variable scalaire $\varphi(\rho)$ avec :

$$
(57) \quad \varphi(\rho) = J(u + \rho \delta_i)
$$

ce qui revient, en fait, à minimiser J(u) par rapport seulement à la
composante u_i du vecteur u. Une itération se compose de la résolution
en séquence des problèmes n°i. L'algorithme s'arrête lorsqu'une
itération entière aura été effectuée sans apporter d'amélioration
notable. On voit donc que la méthode est très, très simple.

Les algorithmes basés sur la méthode séquentielle ne diffèrent que
par la précision avec laquelle les problèmes i sont résolus.

Nous avions proposé en [23] de résoudre chaque problème i très grossièrement : avec des scalaires h_i ayant un sens économique (si u_i est un taux d'intérêt, on prend h_i = 0,25 %) qui sont considérés comme fixés, on se contente de calculer φ(-1), φ(0) et φ(1) et de garder la meilleure valeur. Powell [27] propose de résoudre beaucoup plus précisément chaque problème ; la méthode requiert alors moins d'itérations, mais chaque itération coûte plus cher. Finalement, appliquées à STAR, ces deux méthodes donnent des résultats comparables mais décevants par rapport aux méthodes "de gradient".

Nous avons jugé utile de mentionner cette méthode, bien que
Fair [33] ait déjà montré qu'elle n'est pas compétitive, parce qu'elle
est très simple à programmer soi-même (ce qui peut être utile) et
qu'elle est très économique en place mémoire.

IV.3. Le calcul du gradient réduit par différences finies.

Posons :

$$
(58) \quad h = \max_{1 \le i \le R} |h_i|
$$

et soit $J'_h(u)$ le vecteur de $R^R$ de composantes $J'_{h_i}(u)$ avec :

$$
(59) \quad J'_{h_i}(u) = \frac{1}{h_i} [J(u+\delta_i) - J(u)] , \quad i=1,\dots,R
$$
Si l'on suppose que les fonctions $f_t$ et $j_t$ définissant le modèle et le critère sont continûment différentiables pour tout $u \in U$ et tout $x \in x(U)$, alors on démontre aisément que :

$$ (60) \quad \lim_{h \to 0} J'_h(u) = J'(u) \quad , \quad \forall u \in U $$

Donc, avec des scalaires $h_i$ "petits", en calculant $J'_h(u)$ on obtient une "bonne" approximation du gradient réduit.

Cette technique approchée a l'avantage d'une très grande simplicité mais le calcul (59) nécessite la résolution du modèle pour $R$ valeurs différentes du contrôle ($R = 340$ pour STAR). Nous avons déjà montré que cette méthode, qui est pourtant celle utilisée en général, est bien moins efficace que la méthode de l'état adjoint décrite ci-dessous (cf. [34] et [5]) et nous rappelons ici ce résultat (cf. paragraphe IV. sur les résultats numériques).

### IV.4. Le calcul du gradient réduit par la méthode de l'état-adjoint.

Soit $\delta_t$ l'indicateur défini par :

$$ (61) \quad \delta_t = \begin{cases} 1 & \text{si } t \le T \\ 0 & \text{si } t > T \end{cases} $$

Soient $y_t \in \mathbb{R}^n$ et $w_t \in \mathbb{R}^r$ les vecteurs calculés à partir des gradients en $x$ et $u$ du critère $j(x,u)$ défini en (50) :

$$ (62) \quad y_t = \sum_{k=0}^{p} \delta_{t+k} \frac{\partial j_{t+k}}{\partial x_t} , \quad w_t = \sum_{k=0}^{q} \delta_{t+k} \frac{\partial j_{t+k}}{\partial u_t} , \quad \text{pour } t=1,...,T $$

Enfin, soit $\psi \in \mathbb{R}^N$ la solution du système :

$$ (63) \quad K\psi_t = y_t + \sum_{k=0}^{p} \delta_{t+k} \left( \frac{\partial f_{t+k}}{\partial x_t} \right)^* \psi_{t+k} , \quad \text{pour } t=1,...,T $$

où $A^*$ désigne la matrice transposée de $A$ et où $K$ est la matrice (2). On notera qu'à cause de la notation (61), le système (63) prend pour $t=T$ la forme :

$$ (64) \quad K\psi_T = y_T + \left(\frac{\partial f_T}{\partial x_T}\right)^* \psi_T $$
ce qui permet de calculer $\psi_T$. Pour $t=T-1$, le système (63) prend la forme :

$$ (65) \qquad K\psi_{T-1} = y_{T-1} + \left(\frac{\partial f_{T-1}}{\partial x_{T-1}}\right)^* \psi_{T-1} + \left(\frac{\partial f_T}{\partial x_{T-1}}\right)^* \psi_T $$

ce qui, connaissant $\psi_T$, permet de calculer $\psi_{T-1}$, et ainsi de suite.

Nous appellerons **équation adjointe** le système (63) et nous appellerons **état adjoint** sa solution $\psi$. On vient de voir que la résolution du système adjoint se ramène à la résolution de T systèmes linéaires de dimension n que l'on résout pour $t=T, T-1, T-2, \dots, 2, 1$.

On notera que l'état adjoint dépend à la fois du modèle par l'intermédiaire des matrices $\partial f_{t+k} / \partial x_t$ et du critère par l'intermédiaire du vecteur $y_t$.

**Théorème 6** : supposons que les fonctions $f_t$ et $j_t$ définissant le modèle et le critère soient deux fois continûment différentiables pour tout $u \in U$ et tout $x \in X$. Soit un $u$ quelconque de $U$, $x(u)$ la solution correspondante. Supposons que les vecteurs $y_t$, $w_t$ et les matrices $\partial f_{t+k} / \partial x_t$, $\partial f_{t+k} / \partial u_t$ soient calculées pour le couple $u, x(u)$. Alors le gradient $J'(u)$ du critère (51) est donné par :

$$ (66) \qquad J'(u) = w_t + \sum_{k=0}^{q} \delta_{t+k} \left( \frac{\partial f_{t+k}}{\partial u_t} \right)^* \psi_{t+k}, \quad t=1, \dots, T $$

**Démonstration** :

Soit v un élément quelconque de $\mathbb{R}^R$. Soit h un scalaire suffisamment petit pour que $u+hv \in U$. On pose successivement :

$$ (67) \qquad \begin{gathered} u_h = u+hv \, ; \quad x = x(u) \, ; \quad x_h = x(u_h) \, ; \quad \tilde{x} = (x_h-x)/h \, ; \\ \Delta_h J = [J(u_h) - J(u)] / h. \end{gathered} $$

De la définition (51) du critère J, on déduit :

$$ (68) \qquad \Delta_h J = [j(x_h, u_h) - j(x, u)] / h $$

D'où, d'après la définition (50) du critère j :

$$ (69) \qquad \Delta_h J = \sum_{t=1}^{T} \left[ \sum_{k=0}^{p} \left( \frac{\partial j_t}{\partial x_{t-k}}, \tilde{x}_{t-k} \right)_n + \sum_{k=0}^{q} \left( \frac{\partial j_t}{\partial u_{t-k}}, v_{t-k} \right)_r \right] + O(h) $$
Il est clair que $\tilde{x}_t$ et $v_t$ sont nuls pour tout $t \le 0$. Utilisant successivement les notations (61) puis (62), vérifiera que (69) est équivalent à :

$$ (70) \qquad \Delta_h J = \sum_{t=1}^{T} [(y_t, \tilde{x}_t)_n + (w_t, v_t)_r] + O(h) $$

Tirons $y_t$ du système adjoint (63) et portons dans (70). On obtient :

$$ (71) \qquad \Delta_h J = \sum_{t=1}^{T} (\tilde{x}_t, K\psi_t)_n - \sum_{t=1}^{T} \sum_{k=0}^{p} \delta_{t+k} (\tilde{x}_t, (\frac{\partial f}{\partial x_t})^{*}_{t+k})_{n} + \sum_{t=1}^{T} (w_t, v_t)_r + O(h) $$

Mais $\tilde{x}_t$ étant nul pour tout $t \le 0$, ceci est encore équivalent à :

$$ (72) \qquad \Delta_h J = \sum_{t=1}^{T} \left( K \tilde{x}_t - \sum_{k=0}^{p} \frac{\partial f}{\partial x_{t-k}} \tilde{x}_{t-k}, \psi_t \right)_n + \sum_{t=1}^{T} (w_t, v_t)_r + O(h) $$

Par ailleurs, la dérivation du modèle (1) nous donne, avec la notation (2) :

$$ (73) \qquad K \tilde{x}_t = \sum_{k=0}^{p} \frac{\partial f}{\partial x_{t-k}} \tilde{x}_{t-k} + \sum_{k=0}^{q} \frac{\partial f}{\partial u_{t-k}} v_{t-k} + O(h) $$

Portant (73) dans (72), on obtient :

$$ (74) \qquad \Delta_h J = \sum_{t=1}^{T} \sum_{k=0}^{q} \left( \frac{\partial f}{\partial u_{t-k}} v_{t-k}, \psi_t \right)_n + \sum_{t=1}^{T} (w_t, v_t)_r + O(h) $$

Mais, $v_t$ étant nul pour tout $t \le 0$, ceci équivaut à :

$$ (75) \qquad \Delta_h J = \sum_{t=1}^{T} \left( w_t + \sum_{k=0}^{q} \delta_{t+k} \left( \frac{\partial f}{\partial u_t} \right)^{*}_{t+k} \psi_{t+k} \right)_r + O(h) $$

Comparant à l'expression (66), on voit que :

$$ (76) \qquad \lim_{h \to 0} \frac{1}{h} [J(u_h) - J(u)] = \lim_{h \to 0} \Delta_h J = (J'(u), v)_R $$

Le résultat ayant été obtenu pour un v quelconque, le théorème est démontré.

<figure><img src="image_15.png" /></figure>
Il reste à montrer comment on résout le système adjoint (63).
Nous avons montré que ce système pouvait être résolu pour $t=T$, puis
pour $t=T-1$, etc... Par conséquent, si l'on pose :

$$
(77) \quad \hat{y}_t = y_t + \sum_{k=1}^{p} \delta_{t+k} \left( \frac{\partial f}{\partial x_t} \right)^* \psi_{t+k}, \quad t=1, \dots, T
$$

nous voyons que lors de la résolution de (53) pour la période t, le
vecteur ŷ_t est connu et la résolution se ramène à :

$$
(78) \quad K \psi_t = \hat{y}_t + \left( \frac{\partial f}{\partial x_t} \right)^* \psi_t
$$

La matrice ∂f_t / ∂x_t a la même structure que la matrice d'incidence étudiée au chapitre II (cf. Fig.1, page 11). Par conséquent, pour i = m+s+1,...,n (épilogue) et i = 1,...,l (prologue), la variable ψ_t^i se déduit directement de la connaissance des ψ_t^j, j > i.

On commence donc par résoudre (78) pour i=n, n-1,...,m+s+1 puis
après avoir calculé ψ_t^i pour l < i ≤ m+s, on achève le calcul en
calculant ψ_t^i pour i=l,l-1,...,2,1. Le seul problème est donc le calcul
de ψ_t^i pour i=l+1,...,m+s. Nous poserons :

$$
(79) \quad \psi_t^1 = \{\psi_t^i\}_{i=l+1, \dots, m} ; \quad \psi_t^2 = \{\psi_t^i\}_{i=m+1, \dots, m+s}
$$

Avec les notations introduites au paragraphe III.1, on
vérifiera que les équations adjointes pour i=l+1 à m+s s'écrivent :

$$
(80a) \quad \psi_t^1 = A^* \psi_t^1 + C^* \psi_t^2 + \hat{y}_t^1
$$

$$
(80b) \quad L \psi_t^2 = R^* \psi_t^1 + D^* \psi_t^2 + \hat{y}_t^2
$$

On sait que la matrice A est triangulaire inférieure, donc A*
est triangulaire supérieure. Par conséquent, pour tout vecteur ψ_t^2
donné, le système (80) est récursif. Plus précisément, pour tout ψ_t^2,
on sait calculer ψ_t^1 à l'aide de (80a) et en déduire une erreur, notée
E(ψ_t^2), sur (80b). Il suffit donc de calculer ψ_t^2 tel que E(ψ_t^2) = 0.
D'après (80), il est clair que l'application E est affine et qu'il existe
Λ et μ tels que E(ψ_t^2) = Λψ_t^2 + μ. Le calcul de E(0) donne μ. Calculons
maintenant la matrice Λ. Le système adjoint étant résoluble, il est
clair que la matrice I-A* est régulière, et l'on déduit de (80a) :
$$ (81) \qquad \psi_t^1 = (I - A^*)^{-1} (C^* \psi_t^2 + \hat{y}_t^1) $$

Portant (81) ands (80b), on trouve :

$$ (82) \qquad L\psi_t^2 = B^*(I-A^*)^{-1}(C^*\psi_t^2 + \hat{y}_t^1) + D^*\psi_t^2 + \hat{y}_t^2 $$

D'où le résultat :

$$ (83) \qquad \Lambda = L - B^*(I - A^*)^{-1}C^* - D^* $$

Notant que L est une matrice diagonale, et comparant (83) à (16), on constate que $\Lambda = (\partial\varphi / \partial y)^*$, c'est à dire que la matrice $\Lambda$ est connue : elle est la transposée de la matrice utilisée lors de la résolution du modèle.

La résolution du système adjoint (63) pour une période se ramène donc à la résolution de 2 systèmes récursifs de dimension n et d'un système linéaire de dimension s : le système $\Lambda\psi_t^2 + \mu = 0$.

### IV.5. La méthode de pénalisation.

C'est une méthode très connue et très utilisée en recherche opérationnelle (cf. la bibliographie de [31]) mais qui ne semble pas avoir été utilisée jusqu'à présent en économétrie. L'idée est de considérer les équations du modèle comme des contraintes, les variables endogènes comme des contrôle supplémentaires et d'optimiser un critère augmenté. d'un coût de violation des contraintes. Lorsque le coût marginal de violation des contraintes augmente indéfiniment, la solution du problème pénalisé converge vers la solution du problème initial. L'avantage de la méthode est que le modèle n'a pas besoin d'être résolu à chaque évaluation du critère.

Nous proposons ici, en ne pénalisant que les variables de bouclage, une méthode de pénalisation adaptée à la structure quasi-triangulaire du modèle.

Utilisant la notation (55), on appelle $x^b$ l'ensemble des variables de bouclage et $x^a$ l'ensemble des autres variables endogènes du modèle. Pour des raisons de simplicité théorique, on suppose connu l'ensemble $X^b$ tel que pour tout $u \in U$, les variables de bouclage solution du modèle, $x^b(u)$, sont dans $X^b$. Les $x^b$ étant maintenant des variables de
de contrôle, on appelle v le couple $(x^b, u)$ et V le produit cartésien de $X^b$ et U. Par construction, V est un fermé borné non vide de $R^{S+R}$, avec $S=sT$.

Rappelons que le modèle s'écrit :

$$
(84a) \quad x_t^a = f_t^a(x_t, \dots, x_{t-p}, u_t, \dots, u_{t-q}), \quad t=1, \dots, T
$$

$$
(84b) \quad L x_t^b = f_t^b(x_t, \dots, x_{t-p}, u_t, \dots, u_{t-q}), \quad t=1, \dots, T
$$

Pour tout couple $v = \{x^b, u\}$, (84a) permet de calculer le vecteur $x_t^a$ correspondant par la résolution d'un système récursif de dimension n-s. Par conséquent, (84a) équivaut à la donnée d'une application $v \to \tilde{x}(v)$ qui est **distincte** de l'application $u \to x(u)$ équivalente au modèle (1). En effet, pour un $v$ quelconque, les équations de bouclage ne sont en général pas vérifiées et il y a une erreur notée $g_t(v)$ :

$$
(85) \quad g_t(v) = L x_t^b - f_t^b[\tilde{x}_t(v), \dots, \tilde{x}_{t-p}(v), u_t, \dots, u_{t-q}], \quad t=1, \dots, T
$$

Par construction des applications $u \to x(u)$ équivalente au modèle
et $v \to \tilde{x}(v)$ équivalente à (84a), on a :

$$
(86) \quad [g_t(v) = 0, \forall t] \iff [v = (x^b, u) \text{ et } \tilde{x}(v) = x(u)]
$$

Pour tout $\lambda \in \mathbb{R}^S$, avec $S=sT$, on définit le **lagrangien** du système
par :

$$
(87) \quad L_{\varepsilon}(v, \lambda) = j[\tilde{x}(v), u] + \sum_{t=1}^{T} \sum_{i=1}^{S} \left[ \lambda_{t}^{i} + \frac{1}{\varepsilon} g_{t}^{i}(v) \right] g_{t}^{i}(v)
$$

et le **critère pénalisé**, qui est le cas particulier $\lambda=0$ de (87) :

$$
(88) \quad J_{\varepsilon}(v) = L_{\varepsilon}(v,0) = j[\tilde{x}(v),u] + \frac{1}{\varepsilon} \sum_{t=1}^{T} \sum_{i=1}^{S} [g_t^i(v)]^2
$$

On définit alors deux méthodes de résolution approchée du
problème initial P :

Méthode de pénalisation ordinaire :

On appelle $v_{\varepsilon}$ toute solution du problème pénalisé, c'est à dire tout élément de V vérifiant :

$$
(89) \quad J_{\varepsilon}(v_{\varepsilon}) \le J_{\varepsilon}(v) , \forall v \in V
$$
En s'aidant du fait que V est borné et des techniques développées dans [31], on démontre aisément que pour tout $\epsilon > 0$ petit toute solution de (89) fournit une solution approchée du problème P et plus précisément :

**Théorème 7** : sous les mêmes hypothèses de régularité que pour le problème P, le problème (89) a au moins une solution pour tout $\epsilon > 0$ fixé et la suite $\{v_\epsilon\}$ formée par une quelconque de ses solutions pour chaque $\epsilon$ est l'union de sous-suites $\{v_\epsilon^j\}$ chacune d'entre elles convergeant, lorsque $\epsilon$ tend vers zéro par valeurs positives vers un $\hat{v}_j = \{\hat{x}_j^b, \hat{u}_j\}$ tel que $\hat{u}_j$ est l'une des solutions optimales du problème P et tel que $x(\hat{u}_j)$ soit la solution du modèle (1) associée à $\hat{u}$.

**Méthode du lagrangien augmenté :**

Cette méthode (cf. Rockafellar [35]) est une méthode itérative sur le multiplicateur de Lagrange $\lambda$ de (87). Elle se distingue de la dualité ordinaire par l'addition au lagrangien du terme classique de pénalité, addition rendue indispensable par la non-convexité du critère. A l'itération $k$ de la méthode, $\lambda^k$ étant connu, on appelle $v(\lambda^k)$ tout élément de $V$ qui réalise :

$$ (90) \qquad L_{\epsilon}[v(\lambda^k), \lambda^k] \le L_{\epsilon}(v, \lambda^k) , \quad \forall v \in V $$

puis on passe à l'itération suivante en modifiant $\lambda^k$ par la formule :

$$ (91) \qquad \lambda_t^{k+1} = \lambda_t^k + \frac{2}{\epsilon} g_t[v(\lambda^k)] , \quad t=1, \dots, T $$

Nous admettrons (cf. [35]) qu'il existe un $\epsilon_0$ tel que si $\epsilon$ est choisi fixé dans l'intervalle $[0, \epsilon_0]$, alors le processus (90)-(91) converge.

L'avantage de la méthode de Rockafellar est que $\epsilon$ est fixé alors que dans la pénalisation ordinaire il faut le faire tendre vers zéro. Malheureusement, il faut dès le départ le choisir assez petit ce qui pose des problèmes numériques que nous discuterons au paragraphe IV.

Le problème pour appliquer l'une ou l'autre de ces deux méthodes est le calcul du gradient $J'_\epsilon(v)$ du critère (88) ou du gradient en v du lagrangien (87). Le gradient $J'_\epsilon(v)$ se déduisant de $\partial L_\epsilon(v, \lambda) / \partial v$ en faisant $\lambda=0$ dans toutes les expressions permettant de le calculer, nous nous contenterons d'indiquer le calcul du gradient du lagrangien.
**Théorème 8** : supposons que toutes les fonctions définissant le modèle et le critère soient deux fois continûment différentiables, au moins pour tout $u \in U$ et tout $x \in X$. Soit un $\lambda$ quelconque et un $v \in V$ quelconque, avec $v = \{x^b, u\}$. Posons successivement :

$$ (92) \qquad y_t = \sum_{k=0}^{p} \delta_{t+k} \frac{\partial j_{t+k}}{\partial x_t}, \quad w_t = \sum_{k=0}^{q} \delta_{t+k} \frac{\partial j_{t+k}}{\partial u_t}, \quad t=1, \dots, T $$

$$ (93a) \qquad \psi_t^b = \lambda_t + \frac{2}{\epsilon} g_t(v) $$

$$ (93b) \qquad \psi_t^a = y_t^a + \sum_{k=0}^{p} \delta_{t+k} \left( \frac{\partial f}{\partial x_t^a} \right)^* \psi_{t+k} $$

$$ (94a) \qquad \frac{\partial L}{\partial x_t^b} = y_t^b + \sum_{k=0}^{p} \delta_{t+k} \left( \frac{\partial f}{\partial x_t^b} \right)^* \psi_{t+k} $$

$$ (94b) \qquad \frac{\partial L}{\partial u_t} = w_t + \sum_{k=0}^{q} \delta_{t+k} \left( \frac{\partial f}{\partial u_t} \right)^* \psi_{t+k} $$

Alors, si toutes les expressions ci-dessus sont évaluées pour $v$, le gradient en $v$ du lagrangien (87) est, pour chaque période, égal au couple $\{\frac{\partial L_\epsilon}{\partial x_t^b}, \frac{\partial L_\epsilon}{\partial u_t}\}$ défini en (94).

La démonstration de ce théorème est strictement calquée sur celle du théorème 6.

Il faut noter que contrairement à ce qui se passait pour le calcul du gradient réduit évoqué au paragraphe précédent, le système adjoint (93) est totalement récursif. En effet, on commence par appliquer la formule (93a) donnant $\psi_t^b$ puis on résout le système linéaire (93b) qui est récursif (de façon rétrograde), puisque, rappelons-le, la matrice $(\frac{\partial f^a}{\partial x_t^a})^*$ est triangulaire supérieure.

Pour tout $v = \{x^b, u\}$, le critère se ramène donc à la résolution du système (84a) qui est récursif en $x_t^a$ et le calcul du gradient se ramène à l'évaluation de (93a) puis à la résolution du système linéaire (93b) récursif en $\psi_t^a$.

On conçoit donc, et nous le confirmerons ci-dessous, que chaque itération de la méthode d'optimisation coûtera beaucoup moins de temps de calcul avec l'une des deux méthodes proposées ici qu'avec la méthode du gradient réduit proposée au paragraphe précédent. Reste à savoir,
bien évidemment, si l'augmentation probable (*) du nombre d'itérations nécessaires pour atteindre une précision suffisante va ou non contrebalancer cet avantage.

Les résultats obtenus sur le modèle STAR montrent que la méthode du gradient réduit est finalement plus efficace, mais pas de telle façon que la pénalisation soit à rejeter, compte tenu du fait que l'algorithme est encore grandement améliorable.

## IV.6: Comparaisons numériques.

Nous avons effectué nos essais numériques sur le modèle STAR avec comme critère :

$$ (95) \qquad j(x,u) = \sum_{t=1}^{T} \left[ 5 \cdot 10^{-7} \sum_{i=1}^{r} (u_t^i)^2 + 10^{-4} (BACO_t)^2 \right] $$

où BACO est la balance commerciale que l'on cherche à équilibrer. Nous partons systématiquement de $u=0$ comme approximation initiale, ce qui donne un critère de 85 224. La valeur optimale du critère étant légèrement supérieure à 40, nous estimons que le problème est résolu avec une précision suffisante lorsque le critère est inférieur à 50.

Optimisant le modèle sur 5 périodes l'une ou l'autre des méthodes séquentielles (elles donnent des résultats comparables) proposées en IV.2 permettent d'atteindre un critère de 250 en deux minutes d'un IBM 370/168, ce qui est décevant par rapport aux autres méthodes.

Nous avons testé les deux méthodes de calcul du gradient réduit proposées en IV.3 et IV.4 : la méthode des différences finies (M1) et la méthode de l'état adjoint (M2) dans laquelle les matrices $\partial f_{t+k}^i / \partial x_t^j$ et $\partial f_{t+k}^i / \partial u_t^j$ sont remplacées par des différences finies $[f_{t+k}^i(x_t^j+h) - f_{t+k}^i(x_t^j)] / h$ et $[f_{t+k}^i(u_t^j+h) - f_{t+k}^i(u_t^j)] / h$. Il n'est pas raisonnable d'exiger de l'utilisateur le calcul analytique des expressions des jacobiens, ce qui justifie la méthode employée. Précisons cependant que pour le modèle STAR, l'utilisation du dérivateur formel FORMAC qui génère le code FORTRAN des expressions analytiques du jacobien permet de diviser le temps de calcul de M2 environ par un facteur 10 ce qui rend la méthode de l'état adjoint encore plus efficace.

<footnote>(*) le nombre de variables de contrôle ayant augmenté par rapport à la méthode du gradient réduit.</footnote>
La table ci-dessous donne le temps CPU en FORTRAN H sur un IBM 370/168 nécessaire pour calculer le gradient :

<table><thead><tr><th>T</th><th>M1</th><th>M2</th><th>M1 / M2</th></tr></thead><tbody><tr><td>5</td><td>0,75 sec.</td><td>0,12 sec.</td><td>6,2</td></tr><tr><td>10</td><td>2,95 sec.</td><td>0,24 sec.</td><td>12,3</td></tr><tr><td>20</td><td>11,50 sec.</td><td>0,46 sec.</td><td>25</td></tr><tr><td>30</td><td>26,25 sec.</td><td>0,70 sec.</td><td>37,5</td></tr></tbody></table>

- Table 6 -

T : nombre de périodes de simulation.

M1 : calcul du gradient par différences finies

M2 : calcul du gradient par la méthode de l'état-adjoint, les jacobiens étant calculés par différences finies.

Par conséquent, pour $T=30$ par exemple, la méthode de l'état-adjoint est 37 fois plus rapide que la méthode des différences finies et 370 fois plus rapide si l'on dispose de FORMAC (*). On notera que le temps de calcul de M1 est proportionnel à $T^2$ et le temps de calcul de M2 est proportionnel à T.

La méthode des différences finies, à cause de sa grande simpli-cité, est cependant utile pour vérifier les résultats de la méthode de l'état adjoint au moment du test du programme.

Calculant systématiquement le gradient par la méthode de l'état-adjoint, les jacobiens étant calculés par différences finies, nous avons comparé les méthodes d'optimisation :

a) la méthode de Davidon-Fletcher-Powell appliquée au problème (52), i.e. utilisant le gradient réduit : utilisation du programme VA13AD de Harwell [27];

<footnote>(*) ces résultats numériques sont différents de ceux publiés en [5] car nous utilisons ici un programme de résolution de modèle général alors qu'en [5] le programme avait été écrit spécialement pour STAR, ce qui donnait $M1/M2 = 10,6$ pour $T=30$.</footnote>
b) la méthode du gradient conjugué appliquée au problème (52) :
utilisation du programme VA14AD de Harwell [27] ;

c) la méthode de Davidon-Fletcher-Powell (VA13AD) appliquée au problème pénalisé (89) ;

d) la méthode de Davidon-Fletcher-Powell (VA13AD) appliquée au lagrangien augmenté, c'est à dire l'algorithme (90)-(91).

Disons tout de suite que la méthode b) est nettement moins efficace que la méthode a), ce qui est un résultat très connu en général et qui s'est confirmé pour STAR. Certes la méthode du gradient conjugué est plus économique en réservation mémoire. Avec R égal au nombre de variables de contrôle (R=rT, r dimension du contrôle, T nombre de périodes), VA13AD demande un espace de travail de $R(R+13)/2$ mots et VA14AD demande un espace de travail de 6R mots.
Il semble donc que même avec des modèles plus grands que STAR (*), la conservation en mémoire centrale de l'espace de travail de VA13AD ne soit pas un problème et que la méthode du gradient conjugué soit à rejeter.

En ce qui concerne les méthodes c) et d) il faut noter l'absence d'algorithmes itératifs de résolution du modèle d'où bien des ennuis de divergence possible de ces algorithmes sont évités. Par contre se pose un problème nouveau et difficile : le changement d'unités pour les variables de bouclage indispensable pour le bon conditionnement des problèmes de minimisation (89) et (90). Nous ne connaissons pas actuellement d'autres méthodes que celle des essais successifs avec plusieurs valeurs des nouvelles unités. Les tests effectués sur STAR ont conduits à la constatation suivante : avec le changement de variable $\hat{x}^b = \alpha x^b$, il existe un $\hat{\alpha}$"optimal" tel que : si $\alpha < \hat{\alpha}$, VA13AD appliqué à (89) et (90) ne converge pas et si $\alpha$ augmente en restant supérieur à $\hat{\alpha}$, VA13AD converge mais de plus en plus lentement.

Un deuxième problème est le choix de $\epsilon$. Nous n'avons pas réussi à trouver de valeurs de $\epsilon$ assurant un fonctionnement satisfaisant

<footnote>(*) notons que ce qui compte ici n'est pas le nombre d'équations du modèle, mais le nombre de variables de contrôle.</footnote>
de l'algorithme du lagrangien augmenté, algorithme qui, dans l'état actuel de nos connaissances, est donc à rejeter. Pour la pénalisation ordinaire il faut choisir un $ε_o$ initial et une loi de variation de $ε$ au cours de l'algorithme. Nous avons testé sur le modèle STAR un grand nombre de lois et nous proposons la suivante :

$J^*$ étant la valeur optimale du critère pour le problème initial (52) on montre aisément que, pour tout $ε > 0$, toute solution optimale $u_ε$ du problème (89) vérifie $J_ε(u_ε) \le J$. Nous supposons alors connue une approximation grossière $\hat{J}$ de $J^*$ et dès que l'on a $J_ε(u_ε) \le \hat{J}$, on divise $ε$ par 10 et on recommence. Pour STAR, avec $J = 85221$ pour $u=0$ et $J^* = 40$, nous avons supposé connue l'approximation $\hat{J}=50$. Avec une approximation plus grossière, $\hat{J}=100$ par exemple, l'algorithme converge plus vite mais vers une solution du modèle donnant un critère de l'ordre de $\hat{J}$.

La table ci-dessous donne les temps de calcul sur IBM 370/168 nécessaires pour partir d'un critère de 85221, trouver une solution du modèle avec un critère inférieur à 50, le critère optimal étant égal à 40.

<table><caption>- Table 7 -</caption><tr><td>Méthode a) : gradient réduit</td><td>20 secondes</td></tr><tr><td rowspan="4">Méthode c) : pénalisation ordinaire;<br/>ε initial</td><td>convergence vers<br/>une mauvaise valeur<br/>du critère</td></tr><tr><td>= 10<sup>3</sup> ---- 32 secondes</td></tr><tr><td>= 10<sup>4</sup> ---- 59 secondes</td></tr><tr><td>= 10<sup>5</sup> ---- 38 secondes</td></tr><tr><td></td><td>= 10<sup>6</sup> ---- 47 secondes</td></tr></table>
La conclusion de ce paragraphe est que, au moins en application
à STAR, la combinaison de la méthode du gradient réduit, de la méthode
de l'état-adjoint pour calculer ce gradient et de l'algorithme de
Davidon-Fletcher-Powell est de loin la meilleure méthode d'optimisation.
A notre avis, on peut rejeter la méthode séquentielle (cf. IV.2), le
calcul du gradient par différences finies (cf. IV.3) et la méthode du
gradient conjugué. Les méthodes basées sur la pénalisation des
équations de bouclage sont actuellement moins performantes que le
gradient mais leur implémentation concrète est tellement améliorable
que nous ne pouvons les rejeter.

IV.7. *Traitement des jacobiens du modèle.*

L'utilisation des techniques décrites aux paragraphes IV.4 et
IV.5 exige le calcul, le stockage et la multiplication par les vecteurs
$\psi_{t+k}$ des matrices $((\partial f_{t+k}^i / \partial x_t^j))$ avec $i,j=1,...,m$ et $k = 0,...,K$
où m est le nombre d'équations du modèle et K le nombre maximum de
retard sur les variables endogènes. Cela fait $(K+1)m^2$ éléments à
calculer et stocker avec $K=3$, $m=91$, $(K+1)m^2=33124$ pour STAR ceci en
admettant que l'on ne stocke pas ces matrices pour les différentes
valeurs du temps. Il faut y ajouter 3640 éléments pour les matrices
$((\partial f_{t+k}^i / \partial u_t^j))$.

En réalité, sur les 33124 éléments des matrices $\partial f_{t+k} / \partial x_t$,
235 seulement ne sont pas **identiquement nuls** soit une moyenne de 2,6
variable par second membre d'équation et sur ces 235 éléments non
identiquement nuls, 136 sont des constantes connues et seulement 99
éléments doivent être calculés à chaque itération. De plus, sur les
3640 éléments des matrices $\partial f_{t+k} / \partial u_t$, 10 seulement ne sont pas nuls
et ils sont tous constants.

Ne pas tenir compte de cette structure des jacobiens du modèle
entraîne les inconvenients suivants : dans 99 % des cas on calcule
un élément du jacobien dont on sait a priori qu'il est nul, on
stocke 140 fois trop d'éléments et lorsque l'on résout le système
adjoint, 99 % des opérations sont des multiplications et des additions
dont un des opérandes est nul.
Nous avons donc stocké les jacobiens dans des vecteurs FX(K) et FU(K) de dimensions respectives 235 et 10 en conservant des informations obtenues par différences finies ou par analyse syntaxique du modèle permettant de savoir que le $K^{ème}$ élément de FX est en réalité la dérivée de la $i^{ème}$ équation, prise à t+k, par rapport à $x_t^j$ et permettant de savoir si cet élément du jacobien est constant ou variable (on ne calcule que dans le deuxième cas).

Cette technique, dite des "matrices compressées", nous a permis de réduire considérablement la réservation mémoire nécessaire et surtout d'économiser beaucoup de temps machine en évitant un grand nombre d'opérations inutiles.
Références

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

[1] P. NEPOMIASTCHY ; B. OUDET : "Moduléco, système informatique pour la construction et l'utilisation de modèles macroéconomiques", Bulletin de Liaison de la Recherche en Informatique, n°41, janvier 1978.

[2] D. FOUQUET ; J.M. CHARPIN ; H. GUILLAUME ; P.A. MUET ; D. VALLET : "DMS, Modèle de prévision à moyen terme", Economie et Statistiques, N°79, juin 1976.

[3] A. DRUD : "Method for Control of Complex Dynamic Systems", IMSOR, n°27, Lyngby (DK), 1976 ; Communication au NBER Conference on Control and Economics, New Haven (USA), 25-27 mai 1977.

[4] L. LASDON : Communication au NBER Conference on Control and Economics, New Haven, (USA), 25-27 mai 1977.

[5] P. NEPOMIASTCHY ; A. RAVELLI : "Adapted Methods for Solving and Optimizing Quasi-Triangular econometric models", présenté au NBER Conference on Control and Economics, New Haven (USA) 25-27 mai 1977, et Annals of Economic and Social Measurement, 6/5, 1978.

[6] A. COUTIERE : "un modèle du système monétaire français", Statistiques et Etudes financières, INSEE, 1975/77.

[7] J. MISTRAL : "STAR, Modèle de Prévision à court-moyen terme", Séminaire sur les modèles macroéconomiques (CNRS/NSF), Gif Sur Yvette, 22-26 novembre 1976.

[8] D. FOUQUET : "Les multiplicateurs de DMS", publ. INSEE, 24.1.1977.

[9] B. OUDET : Etude de la dynamique déterministe à court terme des modèles macroéconomiques : application au modèle STAR, Thèse de Sciences Economiques, Univ. de Grenoble I, 1976.

[10] M. DELEAU ; P. MALGRANGE : "Méthodes d'analyse des modèles empiriques", Annales de l'INSEE, n°20, 1975.

[11] "Rapport du groupe de travail : étude des propriétés du modèle MINIDMS", Publ. ENSAE-INSFE-CEPREMAP, 25 juillet 1977.

[12] J. DUFOUR ; G. GILLES : "Approche informationnelle de l'analyse structurale des systèmes complexes", à paraître.

[13] A. BENSOUSSAN ; E. HURST ; B. NASLUND : Management Applications of Modern Control Theory, North-Holland, 1974.
[14] M. ATHANS ; E. KUH ; R.S. PINDYCK : "Stabilization policy using an Imperfect Model and Sequential Control", Communication au NBER Conference on Control and Economics, New Haven (USA), 25-27 mai 1977.

[15] D. KENDRICK : "Applications of Control Theory to Macroeconomics", Annals of Economic and Social Measurement, vol. 5, n°2, 1976, pp. 171-190.

[16] K. WALL ; J.H. WESCOTT : "Policies Optimization Studies Using a Simple UK Economy Control Model", Automatica, Vol. 12, n°6, nov. 76, pp. 551-564.

[17] L.F. PAU : "Research on Optimal Control Adapted to Macro- and Microeconomics : a survey", IFAC 78, 7th World Congress (Finlande) à paraître dans Annals of Economic and Social Measurement.

[18] R.C. FAIR : "A model of Macroeconomic Activity", Vol. II, the Empirical Model, Ballinger Publishing Company, 1976.

[19] K. WALL : NBER Cambridge memo, 26 juin 1976.

[20] E. HELLERMAN ; D.C. RARICK : "The Partitionned Preassigned Pivot Procedure (P4)" in Sparse Matrices and their Applications, D.J. Rose and R.A. Willoughby (eds.), Plenum Press, 1972.

[21] D.V. STEWARD : "Partitionning and Tearing Systems of Equations", SIAM Journal of Numerical Analysis, série B, Vol. 2, 1965.

[22] "A macro-model for the Dutch Economy", Working Group Conference of European Planners, La Haye, 24-25 avril 1969.

[23] P.NEPOMIASTCHY ; B. OUDET : "Comparison of Optimization Techniques for a French Monetary Model", Communication au NBER Conference on Control and Economics, Palo Alto (USA), 26-28 mai 1976.

[24] K. BROWN : "A quadratically Convergent Newton like method Based Upon Gaussian Elimination", SIAM J. on Numerical Analysis, 6, pp. 560-569, 1969.

[25] M. RACHDI : "Thèse de 3e Cycle", à paraître.

[26] B. DEMIDOVITCH ; I. MARON : Éléments de calcul numérique, édit. Mir, trad. française, Moscou, 1973.

[27] Bibliothèque HARWELL, A.E.R.E. Report. Didcot, Oxfordshire (GB).

[28] J.M. ORTEGA ; W.C. RHEINBOLDT : Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970.
[29] C.G. BROYDEN : "A class of Methods for Solving Nonlinear Simultaneous Equations", *Math. Comp.*, 19, pp. 577-593, 1965.

[30] D. GABAY ; D. LUENBERGER : "Efficiently Converging Minimization Methods based on the Reduced Gradient", *Publ. Univ. de Stanford*, 1976.

[31] P. NEPOMIASTCHY : "Méthode de pénalisation et applications", Doc. Travail IERESM (Bruxelles), n°75 - 26 et N° 75-28, 1975.

[32] J. CEA : Optimisation. *Théorie et algorithmes*, Dunod, 1972.

[33] R.C. FAIR : "On the Solution of Optimal Problems as Maximization Problems", *Annals of Economic and Social Measurement*, vol. 3, n°1, Jan. 1974, pp. 135-154.

[34] P. NEPOMIASTCHY : "Méthodes d'optimisation adaptées aux modèles macroéconomiques" *Colloque International sur l'Analyse des Systèmes et ses Orientations Nouvelles*, Déc. 1976, Proceedings publiés par Springer-Verlag, à paraître.

[35] R.T. ROCKAFFELLAR : "Augmented Lagrange Multiplier Functions and Duality in Non Convex Programming", *SIAM J. on Control*, vol. 12, n°2, pp. 268-285, 1974.