Estimation du niveau et de la variation du chômage au moyen de modèles de séries chronologiques structurels
Section 4. Estimation sur petits domaines de séries chronologiques

• Précédent
• Suivant

Les estimations mensuelles initiales du domaine pour les vagues distinctes, accompagnées d’estimations de la variance et de la covariance, sont les données d’entrée pour les modèles de séries chronologiques. À l’étape suivante, on appliquera les MSCS pour lisser les estimations initiales et corriger le biais de renouvellement. On utilise les modèles estimés pour prévoir des fractions du chômage provincial, des tendances du chômage provincial et des variations de ces tendances d’un mois à l’autre. Dans la sous-section 4.1, les MSCS sont définis, puis exprimés en tant que modèles espace-état ajustés dans un cadre fréquentiste. La sous-section 4.2 explique comment ces MSCS peuvent être exprimés sous la forme de modèles multiniveaux de séries chronologiques ajustés dans un cadre bayésien hiérarchique.

4.1  Modèle espace-état

Dans la présente section, on élabore un modèle de séries chronologiques structurel pour les données mensuelles à l’échelle provinciale de 12 provinces simultanément afin de tirer parti des données d’échantillonnage temporelles et transversales. Soit ${\widehat{\overline{Y}}}_{it}\mathrm{<mo>=</mo>}{\left({\widehat{\overline{Y}}}_{it1}\mathrm{,}\dots \mathrm{,}{\widehat{\overline{Y}}}_{it5}\right)}^t$ le vecteur à cinq dimensions contenant les estimations par la régression des données d’enquête ${\widehat{\overline{Y}}}_{itp}$ définies par (3.1) dans la période $t$ et le domaine $i.$ On peut modéliser ce vecteur à l’aide du modèle de séries chronologiques structurel suivant (Pfeffermann, 1991; van den Brakel et Krieg, 2009, 2015) :

$${\widehat{\overline{Y}}}_{it}\mathrm{<mo>=</mo>}{\iota }_5{\theta }_{it}+{\lambda }_{it}+e_{it}\mathrm{,}(4.1)$$

où ${\iota }_5\mathrm{<mo>=</mo>}{\left(\mathrm{1,}\mathrm{1,}\mathrm{1,}\mathrm{1,}1\right)}^t,$ ${\theta }_{it}$ est un scalaire désignant le paramètre de population vrai pour la période $t$ dans un domaine $i,$ ${\lambda }_{it}$ est un vecteur à cinq dimensions qui modélise le biais de renouvellement et $e_{it}$ un vecteur à cinq dimensions avec erreurs d’échantillonnage. Le paramètre de population ${\theta }_{it}$ dans (4.1) est modélisé comme suit :

$${\theta }_{it}\mathrm{<mo>=</mo>}{L}_{it}+S_{it}+{ϵ}_{it}\mathrm{,}(4.2)$$

où $L_{it}$ désigne un modèle de tendance stochastique pour saisir la variation de faible fréquence (tendance plus cycle économique), $S_{it}$ une composante saisonnière stochastique pour modéliser les fluctuations mensuelles et ${ϵ}_{it}$ un bruit blanc pour la variation inexpliquée dans ${\theta }_{it}.$ Pour la composante de tendance stochastique, on utilise le modèle dit de tendance lisse, défini par l’ensemble suivant d’équations :

$$L_{it}=L_{it-1}+R_{it-1}\mathrm{,}{R}_{it}=R_{it-1}+{\eta }_{R\mathrm{,}it}\mathrm{,}{\eta }_{R\mathrm{,}it}\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_{Ri}^2\right)\mathrm{.}(4.3)$$

Pour la composante saisonnière stochastique, on utilise la formule trigonométrique, voir Boonstra et van den Brakel (2016) pour plus de détails. Le bruit blanc dans (4.2) est défini comme étant ${ϵ}_{it}\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_{{ϵ}_i}^2\right).$

Le biais de renouvellement entre les séries d’estimations par la régression des données d’enquête est modélisé dans (4.1) avec ${\lambda }_{it}\mathrm{<mo>=</mo>}{\left({\lambda }_{it1}\mathrm{,}{\lambda }_{it2}\mathrm{,}{\lambda }_{it3}\mathrm{,}{\lambda }_{it4}\mathrm{,}{\lambda }_{it5}\right)}^t.$ On identifie le modèle en prenant ${\lambda }_{it1}\mathrm{<mo>=</mo>\; 0.}$ Cela implique que le biais relatif dans les vagues de suivi par rapport à la première vague est estimé et qu’on suppose que les estimations par la régression des données d’enquête de la première vague sont les approximations les plus fiables de ${\theta }_{it},$ voir la justification dans van den Brakel et Krieg (2009). Les autres composantes modélisent la différence systématique entre la vague $p$ par rapport à la première vague et elles sont modélisées comme des marches aléatoires pour permettre des profils dépendant du temps dans le biais de renouvellement.

$${\lambda }_{itp}\mathrm{<mo>=</mo>}{\lambda }_{it-\mathrm{1\; <mo>;</mo>}p}+{\eta }_{\lambda \mathrm{,}itp}\mathrm{,}{\eta }_{\lambda \mathrm{,}itp}\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_{{\lambda }_i}^2\right)\mathrm{,}p\mathrm{<mo>=</mo>\; 2,}\mathrm{3,}\mathrm{4,}5.(4.4)$$

Enfin, on a élaboré un modèle de séries chronologiques pour les erreurs d’enquête. Soit $e_{it}\mathrm{<mo>=</mo>}{\left(e_{it1}\mathrm{,}{e}_{it2}\mathrm{,}{e}_{it3}\mathrm{,}{e}_{it4}\mathrm{,}{e}_{it5}\right)}^t$ le vecteur à cinq dimensions contenant les erreurs d’enquête des cinq vagues. On utilise les estimations de la variance des estimations par la régression des données d’enquête comme information a priori dans le modèle de séries chronologiques pour tenir compte de l’hétéroscédasticité causée par les tailles d’échantillon variables dans le temps au moyen du modèle d’erreur d’enquête suivant :

$$e_{itp}=\sqrt{v\left({\widehat{\overline{Y}}}_{itp}\right)}{\tilde{e}}_{itp}\mathrm{,}(4.5)$$

et $v\left({\widehat{\overline{Y}}}_{itp}\right)$ défini par (3.2). Comme la première vague est observée pour la première fois, il n’y a pas d’autocorrélation avec des échantillons observés auparavant. Pour modéliser l’autocorrélation entre les erreurs d’enquête des vagues de suivi, on dérive des modèles AR appropriés pour ${\tilde{e}}_{itp},$ en appliquant les équations de Yule-Walker aux coefficients de corrélation

$$\frac{n_{i{t}_1{p}_1{t}_2{p}_2}}{\sqrt{n_{i{t}_1{p}_1}{n}_{i{t}_2{p}_2}}}{\widehat{\rho }}_{t_1{p}_1{t}_2{p}_2}\mathrm{,}(4.6)$$

qui sont dérivés des microdonnées décrites à la section 3. À partir de cette analyse, on suppose un modèle AR(1) pour les vagues 2 à 5 où les coefficients d’autocorrélation dépendent de la vague et du mois. Ces considérations donnent le modèle suivant pour les erreurs d’enquête :

$$\begin{array}{ll}{\tilde{e}}_{it1}\hfill & ={\nu }_{it1}\mathrm{,}{\nu }_{it1}\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_{{\nu }_{i1}}^2\right)\mathrm{,}\hfill \\ {\tilde{e}}_{itp}\hfill & ={\varrho }_{it\left(p-1\right)p}{\tilde{e}}_{i\left(t-3\right)\left(p-1\right)}+{\nu }_{itp}\mathrm{,}{\nu }_{itp}\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_{{\nu }_{ip}}^2\right)\mathrm{,}p\mathrm{<mo>=</mo>\; 2,}\dots \mathrm{,}\mathrm{5,}(4.7)\hfill \end{array}$$

avec ${\varrho }_{it\left(p-1\right)p}$ les coefficients d’autocorrélation partielle dépendant du temps entre les vagues $p$ et $p-1$ dérivés de (4.6). Par conséquent, $\text{Var}\left(e_{it1}\right)\mathrm{<mo>=</mo>}v\left({\widehat{\overline{Y}}}_{it1}\right){\sigma }_{{\nu }_{i1}}^2$ et $\text{Var}\left(e_{itp}\right)\mathrm{<mo>=</mo>}v\left({\widehat{\overline{Y}}}_{itp}\right){\sigma }_{{\nu }_{ip}}^2/\left(1-{\varrho }_{it\left(p-1\right)p}^2\right)$ pour $p\mathrm{<mo>=</mo>\; 2,}\dots \mathrm{,}5.$ Les variances ${\sigma }_{{\nu }_{ip}}^2$ sont des paramètres d’échelle avec des valeurs proches de 1 pour la première vague et proches de $\frac{1}{T}{\displaystyle {\sum }_{t\mathrm{<mo>=</mo>\; 1}}^T}\left(1-{\varrho }_{it\left(p-1\right)p}^2\right)$ pour les autres vagues, où $T$ indique la longueur de la série observée.

Le modèle (4.1) utilise les données d’échantillonnage observées pendant les périodes précédentes dans chaque domaine pour améliorer la précision de l’estimateur par la régression des données d’enquête et il tient compte du biais de renouvellement et de la corrélation sérielle induite par le plan d’échantillonnage à panel rotatif. Pour tirer parti des données d’échantillonnage entre domaines, on peut combiner le modèle (4.1) conçu pour les domaines séparés dans un modèle multivarié unique :

$$\left(\begin{array}{c}{\widehat{\overline{Y}}}_{1t}\\ ⋮\\ {\widehat{\overline{Y}}}_{m_{A}t}\end{array}\right)\mathrm{<mo>=</mo>}\left(\begin{array}{c}{\iota }_5{\theta }_{1t}\\ ⋮\\ {\iota }_5{\theta }_{m_{A}t}\end{array}\right)+\left(\begin{array}{c}{\lambda }_{1t}\\ ⋮\\ {\lambda }_{m_{A}t}\end{array}\right)+\left(\begin{array}{c}{e}_{1t}\\ ⋮\\ e_{m_{A}t}\end{array}\right)\mathrm{,}(4.8)$$

où $m_A$ désigne le nombre de domaines, qui est égal à 12 dans la présente application. Ce cadre multivarié permet d’utiliser l’information de l’échantillon entre les domaines en modélisant la corrélation entre les termes de perturbation des différentes composantes structurelles des séries chronologiques (tendance, saisonnière, biais de renouvellement) ou en définissant les hyperparamètres ou les variables d’état de ces composantes comme étant égaux dans tous les domaines. Dans l’article, les modèles présentant une corrélation transversale entre les termes de perturbation de la pente de la tendance (4.3) sont pris en compte, c’est-à-dire que

$$\text{Cov}\left({\eta }_{R\mathrm{,}it}\mathrm{,}{\eta }_{R\mathrm{,}{i}^{{}^{\prime }}{t}^{{}^{\prime }}}\right)\mathrm{<mo>=</mo>}\{\begin{array}{lll}{\sigma }_{Ri}^2\hfill & \text{si}\hfill & i\mathrm{<mo>=</mo>}{i}^{\prime }\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{et}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}t\mathrm{<mo>=</mo>}{t}^{\prime }\hfill \\ {\varsigma }_{Ri{i}^{{}^{\prime }}}\hfill & \text{si}\hfill & i\ne i^{\prime }\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{et}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}t\mathrm{<mo>=</mo>}{t}^{\prime }\hfill \\ 0\hfill & \text{si}\hfill & t\ne t^{\prime }\hfill \end{array}\mathrm{.}(4.9)$$

La structure de covariance la plus parcimonieuse est une matrice diagonale où tous les domaines partagent une même composante de variance, soit ${\sigma }_{Ri}^2\mathrm{<mo>=</mo>}{\sigma }_R^2$ pour tous les $i$ et ${\varsigma }_{Ri{i}^{{}^{\prime }}}\mathrm{<mo>=</mo>\; 0}$ tous les $i$ et $i^{\prime }.$ Il s’agit de modèles de séries chronologiques structurels dits apparemment non liés, qui constituent une méthode synthétique d’utilisation des données d’échantillonnage entre domaines. Une structure de covariance légèrement plus complexe et réaliste serait une matrice diagonale où chaque domaine a une composante de variance distincte, soit ${\varsigma }_{Ri{i}^{{}^{\prime }}}\mathrm{<mo>=</mo>\; 0}$ pour tous les $i$ et $i^{\prime }.$ Dans ce cas, le modèle emprunte de l’information seulement dans le temps et ne tire pas parti de l’information transversale. La structure de covariance la plus complexe permet une matrice de covariance complète. Une forte corrélation entre les perturbations de pente dans tous les domaines peut produire des tendances cointégrées. Cela implique que les tendances communes $q\mathrm{<mo><</mo>}{m}_A$ sont nécessaires à la modélisation de dynamique des tendances pour les domaines $m_A$ et permettent la spécification des modèles dits à tendance commune (Koopman, Harvey, Doornik et Shephard, 1999; Krieg et van den Brakel, 2012). Les premières analyses de MSCS ont montré que la composante saisonnière et la composante de biais de renouvellement sont indépendantes du temps. Par conséquent, il ne serait pas judicieux de modéliser des corrélations entre les termes de perturbation saisonnière et de biais de renouvellement. Étant donné que les hyperparamètres des paramètres du domaine de population de bruit blanc tendent à être nuls, il semble préférable de supprimer complètement cette composante du modèle, ce qui implique que les corrélations de modélisation entre le bruit de la population ne sont pas prises en compte. Les corrélations entre les erreurs d’enquête pour différents domaines ne sont pas non plus prises en compte, car les domaines sont des régions géographiques dont les échantillons sont tirés indépendamment.

Comme solution de substitution à un modèle avec matrice de covariance complète pour les perturbations de pente, on considère un modèle de tendance qui a un modèle de tendance lisse commun à toutes les provinces plus les composantes de tendance $m_A-1$ qui décrivent l’écart de chaque domaine par rapport à cette tendance globale. Dans ce cas, (4.2) est donné par

$$\begin{array}{ll}{\theta }_{1t}\hfill & \mathrm{<mo>=</mo>}{L}_t+S_{1t}+{ϵ}_{1t}\mathrm{,}\hfill \\ {\theta }_{it}\hfill & \mathrm{<mo>=</mo>}{L}_t+L_{it}^{*}+S_{it}+{ϵ}_{it}\mathrm{,}i\mathrm{<mo>=</mo>\; 2,}\dots \mathrm{,}{m}_A\mathrm{.}(4.10)\hfill \end{array}$$

Ici, $L_t$ est la composante de tendance lisse globale, définie par (4.3), et $L_{it}^{*}$ est l’écart par rapport à la tendance globale pour les domaines distincts, définis comme échelles locales

$$L_{it}^{*}\mathrm{<mo>=</mo>}{L}_{it-1}^{*}+{\eta }_{L\mathrm{,}it}\mathrm{,}{\eta }_{L\mathrm{,}it}\stackrel{\text{ind}}{{\displaystyle \sim }}N\left(\mathrm{0,}{\sigma }_{Li}^2\right)\mathrm{,}(4.11)$$

ou comme des tendances lisses comme dans (4.3). Ces modèles de tendance permettent implicitement des corrélations (positives) entre les tendances des différents domaines.

Les paramètres qu’il faut estimer dans la méthode de modélisation de séries chronologiques sont la tendance et le signal. Ce dernier se définit comme la tendance plus la composante saisonnière. La méthode des séries chronologiques convient particulièrement à l’estimation des variations d’un mois à l’autre. Les profils saisonniers compliquent l’interprétation des variations d’un mois à l’autre des estimations directes et des signaux lissés. C’est pourquoi les variations d’un mois à l’autre sont calculées seulement pour les tendances. En raison de la forte corrélation positive entre les niveaux des périodes consécutives, les erreurs-types des variations d’un mois à l’autre du niveau des tendances sont considérablement plus faibles que celles des variations d’un mois à l’autre des estimations directes, par exemple. La variation d’un mois à l’autre de la tendance est définie comme étant ${\Delta }_{it}\left(1\right)\mathrm{<mo>=</mo>}{L}_{it}-L_{it-1}$ pour les modèles avec tendances séparées pour les domaines ou ${\Delta }_{it}\left(1\right)\mathrm{<mo>=</mo>}{L}_t-L_{t-1}+L_{it}^{*}-L_{it-1}^{*}$ pour les modèles avec une tendance globale et des tendances $m_A-1$ pour l’écart des domaines distincts par rapport à la tendance globale. Cette méthode de modélisation est également utile pour l’estimation des évolutions d’une année à l’autre pour les tendances définies comme étant ${\Delta }_{it}\left(12\right)\mathrm{<mo>=</mo>}{L}_{it}-L_{it-12}$ ou ${\Delta }_{it}\left(12\right)\mathrm{<mo>=</mo>}{L}_t-L_{t-12}+L_{it}^{*}-L_{it-12}^{*}.$ Il est également judicieux de calculer les variations d’une année à l’autre pour les signaux, puisque la plus grande part de la composante saisonnière s’annule. Ces évolutions sont définies d’une manière équivalente aux évolutions de la tendance d’une année à l’autre.

On analyse les modèles de séries chronologiques structurels ci-dessus en les mettant sous une forme dite d’espace-état. On applique ensuite le filtre de Kalman pour ajuster les modèles, où les hyperparamètres inconnus sont remplacés par leurs estimations par le MV. L’analyse est réalisée au moyen d’un logiciel développé dans OxMetrics en combinaison avec les sous-routines de SsfPack 3.0 (Doornik, 2009; Koopman, Shephard et Doornik, 1999, 2008). Les estimations par le MV des hyperparamètres sont obtenues à l’aide de la procédure d’optimisation numérique maxBFGS dans OxMetrics. Voir Boonstra et van den Brakel (2016) pour plus de détails sur la représentation espace-état, l’initialisation du filtre de Kalman et le logiciel utilisé aux fins d’ajustement de ces modèles.

4.2  Modèle multiniveau de séries chronologiques

Pour la description de la représentation des séries chronologiques multiniveau des MSCS, les estimations initiales ${\widehat{\overline{Y}}}_{itp}$ sont combinées dans un vecteur $\widehat{\overline{Y}}\mathrm{<mo>=</mo>}{\left({\widehat{\overline{Y}}}_{111}\mathrm{,}{\widehat{\overline{Y}}}_{112}\mathrm{,}\dots \mathrm{,}{\widehat{\overline{Y}}}_{115}\mathrm{,}{\widehat{\overline{Y}}}_{121}\mathrm{,}\dots \right)}^{\prime },$ c’est-à-dire que l’indice de vague s’exécute plus rapidement que l’indice temporel qui est plus rapide que l’indice de domaine. Le nombre de domaines, de périodes et de vagues est représenté par $m_A,$ $m_T$ et $m_P,$ respectivement. La longueur totale de $\widehat{\overline{Y}}$ est par conséquent $m\mathrm{<mo>=</mo>}{m}_A{m}_T{m}_P\mathrm{<mo>=</mo>}$ 12(domaines)* 72(mois)* 5(vagues) = 4 320. De même, les estimations de la variance $v\left({\widehat{\overline{Y}}}_{itp}\right)$ sont mises dans le même ordre sur la diagonale d’une matrice de covariance $m\times m$ $\Phi .$

La matrice de covariance $\Phi$ n’est pas diagonale en raison des corrélations induites par le plan d’échantillonnage par panel. Il s’agit d’une matrice bande parcimonieuse, et l’ordre du vecteur $\widehat{\overline{Y}}$ est tel qu’il produit une largeur de bande minimale, ce qui est avantageux pour les calculs.

Les modèles multiniveaux envisagés pour la modélisation du vecteur des estimations directes $\widehat{\overline{Y}}$ prennent une forme additive linéaire générale

$$\widehat{\overline{Y}}=X\beta +{\displaystyle \sum _{\alpha }{Z}^{\left(\alpha \right)}{v}^{\left(\alpha \right)}+e,}(4.12)$$

où $X$ est une matrice du plan d’expérience $m\times p$ pour les effets fixes $\beta ,$ et les $Z^{\left(\alpha \right)}$ sont les matrices du plan d’expérience $m\times q^{\left(\alpha \right)}$ pour les vecteurs d’effets aléatoires $v^{\left(\alpha \right)}.$ Ici, la somme sur $\alpha$ s’exécute sur plusieurs termes possibles d’effet aléatoire à différents niveaux, comme une tendance lisse à l’échelle nationale, les tendances d’échelle locale provinciales, le bruit blanc, etc. Cela est expliqué plus en détail ci-dessous. Les erreurs d’échantillonnage $e\mathrm{<mo>=</mo>}{\left(e_{111}\mathrm{,}{e}_{112}\mathrm{,}\dots \mathrm{,}{e}_{115}\mathrm{,}{e}_{121}\mathrm{,}\dots \right)}^{\prime }$ sont considérées comme distribuées normalement comme suit :

$$e\sim N\left(0,\Sigma \right)(4.13)$$

où $\Sigma \mathrm{<mo>=</mo>}{\oplus }_{i\mathrm{<mo>=</mo>\; 1}}^{m_A}{\lambda }_i{\Phi }_i$ avec ${\Phi }_i$ comme matrice de covariance pour les estimations initiales de la province $i,$ et ${\lambda }_i$ un paramètre d’échelle de la variance propre à la province qu’il faut estimer. Comme on l’a décrit à la section 3, les variances par rapport au plan dans $\Phi \mathrm{<mo>=</mo>}{\oplus }_i{\Phi }_i$ sont regroupées sur toutes les provinces et, en raison de la nature discrète des données sur le chômage, elles perdent ainsi une partie de leur dépendance à l’égard du niveau de chômage. On a constaté que l’intégration des facteurs d’échelle de la variance ${\lambda }_i$ permet au modèle de rééchelonner les variances par rapport au plan estimées à un niveau améliorant l’ajustement des données.

Pour décrire le modèle général de chaque vecteur $v^{\left(\alpha \right)}$ d’effets aléatoires, nous supprimons l’exposant $\alpha .$ Chaque vecteur $v$ a $q\mathrm{<mo>=</mo>}dl$ composantes correspondant à $d$ effets autorisés à varier sur $l$ niveaux d’une variable de facteur. En particulier,

$$v\sim N\left(0,A\otimes V\right),(4.14)$$

où $V$ et $A$ sont des matrices de covariance $d\times d$ et $l\times l,$ respectivement. Comme à la section 4.1, la matrice de covariance $V$ peut être paramétrée de trois manières différentes. La plupart du temps, il s’agit d’une matrice de covariance non structurée, c’est-à-dire entièrement paramétrée. Les formes plus parcimonieuses sont $V\mathrm{<mo>=</mo>}\text{diag}\left({\sigma }_{v\mathrm{<mo>;</mo>\; 1}}^2\mathrm{,}\dots \mathrm{,}{\sigma }_{v\mathrm{<mo>;</mo>}d}^2\right)$ ou $V\mathrm{<mo>=</mo>}{\sigma }_v^2{I}_d.$ Si $d\mathrm{<mo>=</mo>\; 1},$ les trois paramétrisations sont équivalentes. La matrice de covariance $A$ décrit la structure de covariance entre les niveaux de la variable de facteur et elle est supposée connue. Il est généralement plus pratique d’utiliser la matrice de précision $Q_A\mathrm{<mo>=</mo>}{A}^{-1},$ car elle est parcimonieuse pour de nombreuses structures courantes de corrélation temporelle et spatiale (Rue et Held, 2005).

4.2.1  Relations entre les représentations espace-état et les représentations multiniveaux des séries chronologiques

On peut représenter une tendance lisse unique comme une ordonnée à l’origine aléatoire $\left(d\mathrm{<mo>=</mo>\; 1}\right)$ variant dans le temps $\left(l\mathrm{<mo>=</mo>}{m}_T\right),$ avec une corrélation temporelle déterminée par une matrice de précision parcimonieuse $Q_A$ de bande $m_T\times m_T$ associée à une marche aléatoire du second ordre (Rue et Held, 2005). Dans ce cas $V\mathrm{<mo>=</mo>}{\sigma }_v^2$ et la matrice du plan d’expérience $Z$ est la matrice des indicateurs $m\times m_T$ pour le mois, c’est-à-dire la matrice avec un seul 1 dans chaque rang pour le mois correspondant et des 0 ailleurs. La parcimonie de $Q_A$ et $Z$ peut être exploitée dans les calculs. La matrice de précision pour la composante de tendance lisse a deux vecteurs singuliers, ${\iota }_{m_T}\mathrm{<mo>=</mo>}\left(\mathrm{1,}\mathrm{1,}\dots \mathrm{,}1\right)$ et ${\left(\mathrm{1,}\mathrm{2,}\dots \mathrm{,}{m}_T\right)}^{\prime }.$ Cela signifie que la spécification correspondante (4.14) est complètement non informative sur le niveau global et la tendance linéaire. Afin d’éviter la non-identifiabilité entre les différents termes du modèle, on peut supprimer le niveau global et la tendance de $v$ en imposant les contraintes $Rv\mathrm{<mo>=</mo>\; 0},$ où $R$ est la matrice $2\times m_T$ avec les deux vecteurs singuliers comme rangs. Le niveau global et la tendance sont ensuite inclus dans le vecteur $\beta$ des effets fixes. Dans la représentation espace-état, on obtient ce modèle en définissant un modèle de tendance (4.3) pour tous les domaines, c’est-à-dire $L_{it}\mathrm{<mo>=</mo>}{L}_t$ et $R_{it}\mathrm{<mo>=</mo>}{R}_t$ pour tous les $i.$ Définir les variables d’état pour la tendance comme étant égales dans les domaines est une méthode très synthétique d’utilisation des données d’échantillonnage provenant d’autres domaines, qui se fonde sur des hypothèses qui, dans la plupart des cas, ne se vérifient pas.

On obtient une tendance lisse pour chaque province avec $d\mathrm{<mo>=</mo>}{m}_A,$ $l\mathrm{<mo>=</mo>}{m}_T,$ et $V$ une matrice de covariance $m_A\times m_A,$ soit diagonale avec un seul paramètre de variance, diagonale avec des paramètres de variance $m_A,$ ou non structurée, c’est-à-dire entièrement paramétrée pour ce qui est des paramètres de variance $m_A$ et des paramètres de corrélation $m_A\left(m_A-1\right)/2.$ Dans ce cas, la matrice du plan d’expérience est $I_{m_A}\otimes I_{m_T}\otimes {\iota }_{m_P}.$ Dans la représentation espace-état, on obtient ces modèles avec le modèle de tendance (4.3) et la structure de covariance (4.9).

Un autre modèle de tendance consisterait en une seule tendance globale lisse (marche aléatoire de second ordre) complétée par une tendance d’échelle locale, c’est-à-dire une marche aléatoire ordinaire (de premier ordre) pour chaque province. Ce dernier peut être modélisé de la manière décrite dans le paragraphe précédent, mais avec une matrice de précision associée à une marche aléatoire de premier ordre. Ce modèle de tendance correspond aux modèles (4.10) et (4.11) dans le contexte espace-état. À la différence de la méthode espace-état, il n’est pas nécessaire de supprimer du modèle une des tendances de la marche aléatoire provinciale aux fins d’identifiabilité. Cela s’explique par le fait que les contraintes de la méthode multiniveau sont imposées de façon à ce que la tendance lisse globale ainsi que toutes les tendances de marche aléatoire provinciale soient à somme nulle dans le temps. Les composantes contraintes correspondent aux ordonnées à l’origine globales et provinciales, qui sont incluses séparément dans le modèle en tant qu’effets fixes, avec l’exclusion d’un effet fixe provincial.

On peut également exprimer les effets saisonniers en termes d’effets aléatoires corrélés (4.14). La composante saisonnière trigonométrique est équivalente au modèle saisonnier de la variable nominale équilibrée (Proietti, 2000; Harvey, 2006), qui correspond à des marches aléatoires de premier ordre dans le temps pour chaque mois, sujettes à une contrainte à somme nulle pour tous les mois. Dans ce cas $d\mathrm{<mo>=</mo>\; 12}$ (saisons), $V\mathrm{<mo>=</mo>}{\sigma }_v^2{I}_{12}$ et $l\mathrm{<mo>=</mo>}{m}_T,$ $Q_A$ étant la matrice de précision d’une marche aléatoire de premier ordre. Les contraintes à somme nulle sur les saisons pour chaque temps, ainsi que les contraintes à somme nulle dans le temps pour chaque marche aléatoire peuvent être imposées comme $Rv\mathrm{<mo>=</mo>\; 0}$ avec $R$ étant la matrice $\left(m_T+12\right)\times 12m_T.$

$$R=\left(\begin{array}{l}{\iota }_{12}^{\prime }\otimes I_{m_T}\hfill \\ I_{12}\otimes {\iota }_{m_T}^{\prime }\hfill \end{array}\right).(4.15)$$

Accompagné d’effets fixes pour chaque saison (encore une fois avec une contrainte à somme nulle imposée), ce terme d’effet aléatoire est équivalent à la saisonnalité trigonométrique. Il peut être étendu à une composante saisonnière pour chaque province, avec un paramètre de variance distinct pour chaque province.

Pour tenir compte du biais de renouvellement, le modèle multiniveau comprend des effets fixes pour les vagues 2 à 5. Facultativement, on peut modéliser ces effets dynamiquement en ajoutant des marches aléatoires dans le temps pour chaque vague. Il faut aussi faire un autre choix, à savoir si les effets fixes et aléatoires sont croisés avec la province.

On peut inclure d’autres effets fixes dans le modèle, par exemple ceux associés aux variables auxiliaires utilisées dans les estimations par la régression des données d’enquête. On pourrait aussi modéliser certaines interactions à effet fixe, par exemple saison  $\times$ province ou vague  $\times$ province, en tant qu’effets aléatoires pour réduire le risque de surajustement.

Enfin, on peut ajouter un terme de bruit blanc au modèle pour tenir compte des variations inexpliquées par domaine et par temps dans le signal.

Le modèle (4.12) peut être considéré comme une généralisation du modèle au niveau du domaine de Fay-Herriot. Le modèle de Fay-Herriot ne comprend qu’un seul vecteur d’effets aléatoires non corrélés sur les niveaux d’une variable à un seul facteur (généralement des régions). Les modèles utilisés dans le présent article comprennent diverses combinaisons d’effets aléatoires non corrélés et corrélés sur des régions et des mois. Des modèles de séries chronologiques multiniveaux étendant le modèle de Fay-Herriot ont été présentés dans Rao et Yu (1994), Datta et coll. (1999), et You (2008). Dans Datta et coll. (1999) et You (2008) des modèles de séries chronologiques sont utilisés avec des effets de domaine indépendants et des marches aléatoires de premier ordre dans le temps pour chaque domaine. Rao et Yu (1994) utilisent un modèle avec effets aléatoires indépendants et un processus autorégressif AR(1) stationnaire au lieu d’un modèle de marche aléatoire dans le temps. You et coll. (2003) ont constaté que le modèle de marche aléatoire correspondait légèrement mieux aux données canadiennes sur le chômage que les modèles AR(1) avec un paramètre d’autocorrélation fixé à 0,5 ou 0,75. Nous ne nous intéresserons pas aux modèles AR(1) dans cet article et nous renvoyons à Diallo (2014) pour une approche permettant à la fois des tendances stationnaires et non stationnaires. Comparé aux références ci-dessus, notre modèle présente une nouvelle caractéristique, à savoir que l’on tient compte des tendances lisses plutôt que des marches aléatoires de premier ordre ou des composantes autorégressives, ou en plus de celles-ci. Nous incluons également les effets aléatoires indépendants de zone dans le temps comme terme de bruit blanc tenant compte de la variation inexpliquée au niveau d’agrégation d’intérêt.

4.2.2  Estimation de modèles multiniveaux de séries chronologiques

On utilise une approche bayésienne pour adapter le modèle (4.12)-(4.14). Cela signifie que nous avons besoin de distributions a priori pour tous les (hyper)paramètres du modèle. Les lois a priori suivantes sont utilisées :

• On attribue aux paramètres de variance au niveau des données ${\lambda }_i$ pour $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{m}_A$ des lois du chi carré inverses avec des degrés de liberté et des paramètres d’échelle égaux à 1.
• Les effets fixes se voient attribuer une loi a priori normale avec une matrice de variance diagonale fixe et une moyenne nulle, présentant de très grandes valeurs (1e10).
• Pour une matrice de covariance entièrement paramétrée $V$ dans (4.14), nous utilisons la loi de de Wishart inverse a priori mise à l’échelle proposée dans O’Malley et Zaslavsky (2008) et recommandée par Gelman et Hill (2007). Conditionnellement à un paramètre vectoriel dimensionnel $d$ $\xi ,$

$$V|\xi \sim \text{Inv}-\text{Wishart}\left(V|v,\text{diag}\left(\xi \right)\Psi \text{diag}\left(\xi \right)\right)(4.16)$$

• où $\nu \mathrm{<mo>=</mo>}d+1$ est choisi et $\Psi \mathrm{<mo>=</mo>}{I}_d.$ On attribue au vecteur $\xi$ une distribution normale $N\left(\mathrm{0,}{I}_d\right).$
• On attribue à tous les autres paramètres de variance apparaissant dans une matrice diagonale $V$ dans (4.14), conditionnellement à un paramètre auxiliaire $\xi ,$ des lois a priori du chi carré inverses avec un degré de liberté de 1 et un paramètre d’échelle ${\xi }^2.$ On attribue à chaque paramètre $\xi$ une loi a priori $N\left(\mathrm{0,}1\right).$ Marginalement, les paramètres de l’écart-type ont des demi-distributions de Cauchy a priori. Gelman (2006) montre que ces lois a priori sont de meilleures lois a priori par défaut que les lois a priori du chi carré inverses, plus courantes.

Le modèle est ajusté au moyen d’un échantillonnage par méthode de Monte-Carlo par chaînes de Markov (MCMC), en particulier l’échantillonneur de Gibbs (Geman et Geman, 1984; Gelfand et Smith, 1990). Les modèles multiniveaux considérés appartiennent à la classe des modèles gaussiens latents additifs dont les termes d’effet aléatoire sont des champs aléatoires gaussiens de Markov (GMRF), et nous utilisons la matrice parcimonieuse et les techniques d’échantillonnage par blocs décrites dans Rue et Held (2005) pour ajuster efficacement ces modèles aux données. De plus, la paramétrisation selon les paramètres auxiliaires susmentionnés $\xi$ (Gelman, Van Dyk, Huang et Boscardin, 2008) améliore considérablement la convergence de l’échantillonneur de Gibbs utilisé. Voir Boonstra et van den Brakel (2016) pour plus de détails sur l’échantillonneur de Gibbs utilisé, y compris les spécifications des distributions conditionnelles complètes. Les méthodes sont mises en œuvre en R au moyen du paquet en R mcmcsae (Boonstra, 2016).

Pour chaque modèle considéré, l’échantillonneur de Gibbs est exécuté dans trois chaînes indépendantes avec des valeurs de départ aléatoires. Chaque chaîne s’exécute pendant 2 500 itérations. Les 500 premiers tirages sont supprimés en tant « qu’échantillon de rodage ». Sur les 2 000 tirages restants de chaque chaîne, nous conservons un tirage sur cinq pour économiser de la mémoire tout en réduisant l’effet d’autocorrélation entre les tirages successifs. Cela nous laisse 3 * 400 = 1 200 tirages pour calculer les estimations et les erreurs-types. On a constaté que le nombre réel de tirages indépendants était de près de 1 200 pour la plupart des paramètres du modèle, ce qui signifie que la plupart des autocorrélations sont éliminées par la réduction. On évalue la convergence de la simulation par MCMC au moyen de tracés et de graphiques d’autocorrélation ainsi que du facteur de réduction d’échelle potentiel de Gelman-Rubin (Gelman et Rubin, 1992), qui diagnostique le mélange de chaînes. Les diagnostics suggèrent que toutes les chaînes convergent correctement à l’étape de rodage, et que les chaînes se mélangent bien, puisque tous les facteurs de Gelman-Rubin sont proches de 1. De plus, les erreurs de simulation de Monte-Carlo estimées (qui tiennent compte de toute autocorrélation restant dans les chaînes) sont faibles comparativement aux erreurs-types a posteriori pour tous les paramètres, si bien que le nombre de tirages conservés suffit à nos besoins.

On peut exprimer les paramètres à estimer d’intérêt sous forme de fonctions des paramètres, et l’application de ces fonctions aux données de sortie de la méthode MCMC des paramètres donne des résultats tirés des lois a posteriori pour ces paramètres à estimer. Dans l’article, nous synthétisons ces tirages en fonction de leur écart moyen et de leur écart-type, qui servent respectivement d’estimations et d’erreurs-types. Toutes les estimations prises en compte peuvent être exprimées sous forme de prédicteurs linéaires, c’est-à-dire sous forme de combinaisons linéaires des paramètres du modèle. On calcule les estimations et les erreurs-types pour les paramètres à estimer suivants :

• Signal : le vecteur ${\theta }_{it},$ y compris tous les effets fixes et aléatoires, sauf ceux associés aux vagues 2 à 5. Ces valeurs correspondent aux valeurs ajustées $X\beta +{\displaystyle {\sum }_{\alpha }{Z}^{\left(\alpha \right)}{v}^{\left(\alpha \right)}}$ associées à une rangée sur cinq $\mathrm{1,}\mathrm{6,}\mathrm{11,}\dots$ de $\widehat{\overline{Y}}$ et des matrices du plan d’expérience.
• Tendance : prédiction de la tendance à long terme. On réalise ce calcul en incorporant seulement les composantes des tendances de chaque modèle dans le prédicteur linéaire. Pour la plupart des modèles considérés, la tendance correspond aux chiffres désaisonnalisés, c’est-à-dire à des prédictions du signal dont tous les effets saisonniers sont éliminés.
• Croissance de la tendance : les différences entre les tendances au cours de deux mois consécutifs.

• Précédent
• Page principale de l'article
• Suivant