Exemple de Régression par les Moindres Carrés Généralisés (GLS)
Exemple de Régression par les Moindres Carrés Généralisés (GLS)#
La méthode canonique du paquet statsmodels pour l’ajustement des modèles de régression linéaire selon le critère des moindres carrés généralisés (GLS) est la méthode GLS.
Les paramètres requis sont :
endog:array-likeUne variable de réponse endogène sous forme d’objet compatible avec les arrays à une dimension
numpy.exog:array-likeUn tableau \(n\times k\) où \(n\) est le nombre d’observations et \(k\) est le nombre de régresseurs. Un terme intercept n’est pas inclus par défaut et doit être spécifié par l’utilisateur (au moyen de
add_constant).sigma:scalarouarrayUn tableau ou un scalaire de type
numpydésignant la matrice de variance-covariance pondérée \(\Sigma\). La valeur par défaut estNone. Sisigmaest un scalaire, il est supposé quesigmaest une matrice diagonale \(n\times n\) avec le scalaire donné comme valeur de chaque élément diagonal. Sisigmaest un vecteur de longueur \(n\), alorssigmaest supposé être une matrice diagonale avec la valeur donnée sur la diagonale.
On commence par charger les paquets nécessaires, principalement matplotlib.pyplot pour le traçage, et numpy pour la manipulation des matrices et vecteurs et la génération aléatoire.
import numpy as np
import statsmodels.api as sm
Note
statsmodels offre la possibilité de charger des jeux de données classiques identiques à ceux du logiciel R, via le sous-module datasets pour réaliser des tests, comparer des modèles étudiés, ou pour des tutoriels d’apprentissage.
Le jeu de données qu’on utilisera pour cette manipulation est Longley. C’est une réalisation de séries temporelles de diverses variables macroéconomiques américaines connues pour être fortement colinéaires. On en extrait la variable de réponse.
data = sm.datasets.longley.load()
data.exog = sm.add_constant(data.exog)
n=len(data.exog)
print(data.exog.head())
const GNPDEFL GNP UNEMP ARMED POP YEAR
0 1.0 83.0 234289.0 2356.0 1590.0 107608.0 1947.0
1 1.0 88.5 259426.0 2325.0 1456.0 108632.0 1948.0
2 1.0 88.2 258054.0 3682.0 1616.0 109773.0 1949.0
3 1.0 89.5 284599.0 3351.0 1650.0 110929.0 1950.0
4 1.0 96.2 328975.0 2099.0 3099.0 112075.0 1951.0
Pour déterminer le paramètre sigma, on estime d’emblée les résidus qui feront figure des erreurs \((\varepsilon_i)_i\), avec la régression par le critère des moindres carrés OLS.
Mise en garde
Cette technique d’estimation de la matrice \(\Sigma\) est dite Feasible Generalized Least Squares (FGLS), et se prête mieux aux échantillons de très grande taille.
ols_resid = sm.OLS(data.endog, data.exog).fit().resid
Selon la documentation du jeu de données Longley, les termes d’erreur suivent un processus stochastiques de type \(\mathrm{AR}(1)\), avec une tendance : \(\varepsilon_i = \beta_0 + \rho\varepsilon_{i-1} + \eta_i\) où \(\eta \sim \mathcal N(\mathbf 0,\Sigma^2)\). \(\rho\) est par conséquent l’autocorrélation des résidus, et peut être facilement estimée en estimant les paramètres de la régression des résidus par les mêmes résidus retardés (\(\varepsilon_i\sim\varepsilon_{i-1}\)).
resid_fit = sm.OLS(np.asarray(ols_resid)[1:], sm.add_constant(np.asarray(ols_resid)[:-1])).fit()
rho = resid_fit.params[1]
print(rho)
-0.3634294908796553
Puisque les termes d’un processus \(\mathrm{AR}(1)\) ont une forte corrélation avec les termes voisins, la matrice de variance-covariance pondérée \(\Sigma\) peut être déterminée comme suit :
sigma = rho**np.array([[abs(j-i) for j in range(n) ]for i in range(n)])
gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)
gls_results = gls_model.fit()
print(gls_results.summary())
GLS Regression Results
==============================================================================
Dep. Variable: TOTEMP R-squared: 0.998
Model: GLS Adj. R-squared: 0.997
Method: Least Squares F-statistic: 724.0
Date: Thu, 05 Jan 2023 Prob (F-statistic): 1.48e-11
Time: 23:24:00 Log-Likelihood: -107.50
No. Observations: 16 AIC: 229.0
Df Residuals: 9 BIC: 234.4
Df Model: 6
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const -3.798e+06 6.71e+05 -5.663 0.000 -5.32e+06 -2.28e+06
GNPDEFL -12.7656 69.431 -0.184 0.858 -169.829 144.298
GNP -0.0380 0.026 -1.448 0.182 -0.097 0.021
UNEMP -2.1869 0.382 -5.719 0.000 -3.052 -1.322
ARMED -1.1518 0.165 -6.970 0.000 -1.526 -0.778
POP -0.0681 0.176 -0.386 0.709 -0.467 0.331
YEAR 1993.9529 342.635 5.819 0.000 1218.860 2769.046
==============================================================================
Omnibus: 1.365 Durbin-Watson: 2.534
Prob(Omnibus): 0.505 Jarque-Bera (JB): 0.885
Skew: 0.209 Prob(JB): 0.642
Kurtosis: 1.926 Cond. No. 5.61e+09
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 5.61e+09. This might indicate that there are
strong multicollinearity or other numerical problems.
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/stats/_stats_py.py:1477: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16
warnings.warn("kurtosistest only valid for n>=20 ... continuing "
La qualité d’ajustement du modèle linéaire par le critère des moindres carrés généralisés, i.e. \(R^2=99\%\), confirme l’origine des données du jeu Longley.