Cet article vise à essayer le modèle de régression à coefficient variable dans le temps, qui est l'une des variantes du modèle d'espace d'états. Comme il est facile d'utiliser R dlm, cette fois je vais appeler et exécuter R dlm de Python via rpy2. Pour plus d'informations sur l'utilisation du modèle de coefficient variable dans le temps dlm avec R, voir Logics of Blue "[modèle de coefficient variable basé sur le temps dlm](http://logics-of-blue.com/dlm%E3%81%AB%E3%82%". 88% E3% 82% 8B% E6% 99% 82% E5% A4% 89% E4% BF% 82% E6% 95% B0% E3% 83% A2% E3% 83% 87% E3% 83% AB / ) »A été très utile. Dans cet article, j'ai changé les variables explicatives de un à deux de ce modèle et j'ai essayé de l'utiliser de Python à rpy2.

1. Qu'est-ce qu'un modèle de régression à coefficient variable dans le temps?

Comme son nom l'indique, il s'agit d'un modèle dans lequel le coefficient de régression change avec le temps.

y_t = a_t x_t + b_t + v_t \quad \cdots (1)

$ y_t $: variable dépendante au temps $ t $ $ x_t $: variable explicative pour le temps $ t $ $ a_t $: coefficient de régression au point temporel $ t $ $ b_t $: Section de temps $ t $ $ v_t $: Erreur d'observation au point temporel $ t $

Comme indiqué, le coefficient de régression change avec le temps. ce temps,

a_t = a_{t-1} + e_t \quad \cdots (2)

Comme indiqué ci-dessus, le coefficient de régression de la période précédente a une erreur de $ e_t $, et le coefficient de régression de la période $ t $ est obtenu. En termes de modèle d'espace d'états, (1) est le modèle d'observation et (2) est le modèle du système.

Cette fois, je voudrais essayer le cas où il y a deux variables explicatives, donc la série d'équations est la suivante.

\begin{aligned}
y_t &= a^{(1)}_t x^{(1)}_t + a^{(2)}_t x^{(2)}_t +b_t + v_t \\
a^{(1)}_t &= a^{(1)}_{t-1} + e^{(1)}_t \\
a^{(2)}_t &= a^{(2)}_{t-1} + e^{(2)}_t \\
b_t &= b_{t-1} + w_t \\
\\
v_t &\sim N(0, \sigma_{v})\\
e^{(1)}_t &\sim N(0, \sigma_{e}^{(1)}) \\
e^{(2)}_t &\sim N(0, \sigma_{e}^{(2)}) \\
w_t &\sim N(0, \sigma_{w})\\
\end{aligned}

Il existe un modèle d'observation et trois modèles de système. Il existe trois variables latentes, $ a ^ {(1)} _t, a ^ {(2)} _t et b_t $ x $ 3t $ d'heures $ t $.

2. Créer des données

Après avoir défini divers paramètres selon l'équation ci-dessus, générez des nombres aléatoires et créez des données. Je me suis référé à [1] pour les spécifications des paramètres de ces données artificielles.

rd.seed(71)

T = 500
v_sd = 20 #Écart type de l'erreur d'observation
i_sd = 10 #écart type de l'interception
a10 = 10
e_sd1 = .5 #Écart type de la variation du coefficient de régression 1
a20 = 20
e_sd2 = .8 #Écart type de la variation du coefficient de régression 1

#Deux coefficients de régression variant dans le temps
e1 = np.random.normal(0, e_sd1, size=T)
a1 = e1.cumsum() + a10
e2 = np.random.normal(0, e_sd2, size=T)
a2 = e2.cumsum() + a20

# intercept
intercept = np.cumsum(np.random.normal(0, i_sd, size=T))

#Variable explicative
x1 = rd.normal(10, 10, size=T)
x2 = rd.normal(10, 10, size=T)

#Variable dépendante
v = np.random.normal(0, v_sd, size=T) #Erreur d'observation
y = intercept + a1*x1 + a2*x2 + v
y_noerror = intercept + a1*x1 + a2*x2  #Y lorsqu'il n'y a pas eu d'erreur d'observation

Le graphique résultant est ci-dessous. Tout d'abord, un graphique des trois variables latentes que vous souhaitez estimer.

Graphique du coefficient de régression 1
Graphique du coefficient de régression partie 2
Graphique en coupe

Les variables explicatives sont des données qui suivent deux distributions normales indépendantes.

Variables explicatives $ x ^ {(1)} _ t, x ^ {(2)} _ t $
Résultat $ y $

En conséquence, le graphique de $ y $ ne semble pas avoir de structure avec juste des nombres aléatoires à l'œil humain. Je voudrais confirmer si les variables latentes peuvent être bien estimées à partir de ces données.

3. Estimation du modèle à l'aide de dlm

Eh bien, voici la production. Je voudrais appeler dlm de R de Python à rpy2 pour estimer le modèle de régression à coefficient variant dans le temps.

3-1. Importation rpy2

Importez rpy2.

import rpy2
from rpy2.robjects.packages import importr
import rpy2.robjects as robjects
from rpy2.rinterface import R_VERSION_BUILD
from rpy2.robjects import pandas2ri
pandas2ri.activate()

Importez la bibliothèque R, dlm.

dlm = importr("dlm")

3-2. Spécification du modèle du modèle d'espace d'états

Ici, vous spécifiez le type de modèle d'espace d'états. Cette fois, puisqu'il s'agit d'un modèle de régression avec des variables explicatives, nous utiliserons dlmModReg. «X» est une variable explicative, et cette fois il y en a deux, «x1» et «x2», alors définissez ces données. dV est la variance de l'erreur du modèle d'observation et représente la variance de $ v_t $ $ \ sigma_v $. «dW» est la distribution des erreurs de modèle de système. Cette fois, il existe trois modèles système liés à $ e ^ {(1)} _t, e ^ {(2)} _t, w_t $.

Cela correspond respectivement à $ \ sigma_ {e} ^ {(1)}, \ sigma_ {e} ^ {(2)}, \ sigma_w $.

`python`


buildDlmReg = lambda theta: dlm.dlmModReg(
    X=df[["x1", "x2"]], 
    dV=np.exp(theta[0]), 
    dW=[np.exp(theta[1]), np.exp(theta[2]), np.exp(theta[3])]
  )

Maintenant que le type de modèle a été décidé, définissez la valeur observée «y» et trouvez les paramètres par la méthode la plus probable. Puisque l'estimation peut être incorrecte en fonction de la valeur initiale, elle est estimée en utilisant dlmMLE deux fois. Cela signifie que le résultat du premier «dlmMLE» est utilisé comme valeur initiale et que le second «dlmMLE» est exécuté pour améliorer la stabilité.

3-3. Estimation de la variance de l'erreur par la méthode la plus probable

`python`


#Effectuer l'estimation la plus probable en deux étapes
%time parm = dlm.dlmMLE(y, parm=[2, 1, 1, 1], build=buildDlmReg, method="L-BFGS-B") 
par = np.array(get_method(parm, "par"))
print("par:", par)
%time fitDlmReg = dlm.dlmMLE(y, parm=par, build=buildDlmReg, method="SANN")

Affiche le résultat.

`python`


show_result(fitDlmReg)

Le par suivant correspond aux résultats de l'estimation des quatre écarts types de $ \ sigma_v, e ^ {(1)} _t, e ^ {(2)} _t, w_t $, mais il est défini sur dlmModReg. Parfois, je mets np.exp entre les deux pour éviter que la valeur ne devienne négative, donc si je prends le résultat np.exp et que je prends ensuite np.sqrt, ce sera l'écart type.

`out`


par [1]  5.9998846  4.5724213 -1.6572242 -0.9232988

value [1] 2017.984

counts function gradient 
   10000       NA 

convergence [1] 0

message NULL

Lorsque vous calculez réellement, vous pouvez voir que les paramètres spécifiés lors de la création des données artificielles sont presque reproduits.

`python`


par = np.asanyarray(get_method(fitDlmReg, "par"))
modDlmReg = buildDlmReg(par)
estimated_sd = np.sqrt(np.exp(np.array(get_method(fitDlmReg, "par"))))
pd.DataFrame({"estimated_sd":estimated_sd, "sd":[v_sd, i_sd, e_sd1, e_sd2]})

	estimated_sd	sd
0	20.084378	20.0
1	9.837589	10.0
2	0.436655	0.5
3	0.630243	0.8

3-4. Exécution du filtrage et du lissage

À présent, en utilisant ce résultat, recherchez la valeur filtrée par le filtre de Kalman, puis recherchez la valeur lissée à partir de cette valeur de filtre.

`python`


#Filtre de Kalman
filterDlmReg = dlm.dlmFilter(y, modDlmReg)

#Lissage
smoothDlmReg = dlm.dlmSmooth(filterDlmReg)

Les résultats de filtrage et de lissage obtenus sont renvoyés dans le monde Python et stockés en tant que ndarray.

`python`


#Obtenir des valeurs filtrées
filtered_a = np.array(get_method(filterDlmReg, "a"))

#Obtenez une valeur lissée
smoothed_a = np.array(get_method(smoothDlmReg, "s"))

3-5. Visualisation des résultats des estimations

Ci-dessous, un graphique des deux résultats, section ʻinterception`, coefficient de régression $ a ^ {(1)} _ t, a ^ {(2)} _ t $. La valeur d'origine des données artificielles est affichée en bleu, la valeur filtrée est affichée en vert et la valeur lissée est affichée en rouge, mais je pense que le mouvement des données d'origine peut être reproduit considérablement. Bien sûr, dlm ne reçoit que trois données, «x1, x2, y», et le type du modèle d'espace d'état, mais la section et le coefficient de régression peuvent être restaurés correctement! (La valeur est rampante seulement près de la valeur initiale)

`python`


plt.figure(figsize=(12,5))
plt.plot(df.intercept, label="intercept(original)")
plt.plot(filtered_a[1:,0], label="intercept(filtered)")
plt.plot(smoothed_a[1:,0], label="intercept(smoothed)")
plt.legend(loc="best")
plt.title("plot of intercept")

`python`


plt.figure(figsize=(12,5))
plt.plot(df.a1, label="a1(original)")
plt.plot(filtered_a[1:,1], label="a1(filtered)")
plt.plot(smoothed_a[1:,1], label="a1(smoothed)")
plt.legend(loc="best")
plt.title("plot of regression coeff a1")

`python`


plt.figure(figsize=(12,5))
plt.plot(df.a2, label="a2(original)")
plt.plot(filtered_a[1:,2], label="a2(filtered)")
plt.plot(smoothed_a[1:,2], label="a2(smoothed)")
plt.legend(loc="best")
plt.title("plot of regression coeff a2")

Enfin, comparons celui sans l'erreur d'observation de la valeur y et celui calculé avec y en utilisant la variable latente lissée. Vous pouvez voir que la vraie valeur y peut être reproduite considérablement en supprimant ici également l'erreur d'observation!

`python`


estimatedLevel = smoothed_a[1:,0] + df.x1*smoothed_a[1:,1] + df.x2*smoothed_a[1:,2]

`python`


plt.figure(figsize=(12,5))
plt.plot(y_noerror, lw=4, label="y(original)")
plt.plot(estimatedLevel, "r", lw=1, label="y(estimated) [no observation error]")
plt.legend(loc="best")
plt.title("plot of target value.")

référence

[1] Logiques du bleu "Modèle à coefficient variant dans le temps par dlm" http://logics-of-blue.com/dlm%E3%81%AB%E3%82%88%E3%82%8B%E6%99%82%E5%A4%89%E4%BF%82%E6%95%B0%E3%83%A2%E3%83%87%E3%83%AB/

[2] Ce code est téléchargé (github) https://github.com/matsuken92/Qiita_Contents/blob/master/General/dlm_time_varying_regression.ipynb

[2] Version de ce code qui peut être supportée par n variables (github) https://github.com/matsuken92/Qiita_Contents/blob/master/General/dlm_time_varying_regression_n-var.ipynb

Appelez dlm depuis python pour exécuter un modèle de régression à coefficient variable dans le temps

1. Qu'est-ce qu'un modèle de régression à coefficient variable dans le temps?

2. Créer des données

3. Estimation du modèle à l'aide de dlm

3-1. Importation rpy2

3-2. Spécification du modèle du modèle d'espace d'états

python

3-3. Estimation de la variance de l'erreur par la méthode la plus probable

python

python

out

python

3-4. Exécution du filtrage et du lissage

python

python

3-5. Visualisation des résultats des estimations

python

python

python

python

python

référence

`python`

`python`

`python`

`out`

`python`

`python`

`python`

`python`

`python`

`python`

`python`

`python`