Rufen Sie dlm von Python aus auf, um ein zeitvariables Koeffizientenregressionsmodell auszuführen

Dies ist ein Artikel zum Testen des zeitvariablen Koeffizientenregressionsmodells, das eine der Variationen des Zustandsraummodells ist. Da es einfach ist, R dlm zu verwenden, werde ich dieses Mal R dlm von Python über rpy2 aufrufen und ausführen. Informationen zur Verwendung des zeitvariablen Koeffizientenmodells von dlm nach R finden Sie unter Logik von Blau "[Zeitvariables Koeffizientenmodell nach 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 / ) ”War sehr hilfreich. In diesem Artikel habe ich die erklärenden Variablen dieses Modells von eins auf zwei geändert und versucht, sie von Python über rpy2 zu verwenden.

1. Was ist ein zeitvariables Koeffizientenregressionsmodell?

Wie der Name schon sagt, handelt es sich um ein Modell, bei dem sich der Regressionskoeffizient mit der Zeit ändert.

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

$ y_t $: Abhängige Variable zum Zeitpunkt $ t $ $ x_t $: Erklärende Variable für die Zeit $ t $ $ a_t $: Regressionskoeffizient zum Zeitpunkt $ t $ $ b_t $: Zeitabschnitt $ t $ $ v_t $: Beobachtungsfehler zum Zeitpunkt $ t $

Wie gezeigt, ändert sich der Regressionskoeffizient mit der Zeit. diese Zeit,

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

Wie oben gezeigt, hat der Regressionskoeffizient der vorherigen Periode einen Fehler von $ e_t $, und der Regressionskoeffizient der Periode $ t $ wird erhalten. In Bezug auf das Zustandsraummodell ist (1) das Beobachtungsmodell und (2) das Systemmodell.

Dieses Mal möchte ich den Fall ausprobieren, in dem es zwei erklärende Variablen gibt, daher ist die Reihe der Gleichungen wie folgt.

\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}

Es gibt ein Beobachtungsmodell und drei Systemmodelle. Es gibt drei latente Variablen: $ a ^ {(1)} _t, a ^ {(2)} _t und b_t $ x $ 3t $ Stunden $ t $.

2. Daten erstellen

Nachdem Sie verschiedene Parameter gemäß der obigen Gleichung eingestellt haben, generieren Sie Zufallszahlen und erstellen Sie Daten. Ich habe auf [1] für die Parameterspezifikationen dieser künstlichen Daten verwiesen.

rd.seed(71)

T = 500
v_sd = 20 #Standardabweichung des Beobachtungsfehlers
i_sd = 10 #Standardabweichung des Abschnitts
a10 = 10
e_sd1 = .5 #Standardabweichung der Variation des Regressionskoeffizienten 1
a20 = 20
e_sd2 = .8 #Standardabweichung der Variation des Regressionskoeffizienten 1

#Zwei zeitvariable Regressionskoeffizienten
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))

#Erklärende Variable
x1 = rd.normal(10, 10, size=T)
x2 = rd.normal(10, 10, size=T)

#Abhängige Variable
v = np.random.normal(0, v_sd, size=T) #Beobachtungsfehler
y = intercept + a1*x1 + a2*x2 + v
y_noerror = intercept + a1*x1 + a2*x2  #Y, wenn kein Beobachtungsfehler aufgetreten ist

Das resultierende Diagramm ist unten. Zunächst ein Diagramm der drei latenten Variablen, die Sie schätzen möchten.

• Diagramm des Regressionskoeffizienten 1

• Diagramm des Regressionskoeffizienten Teil 2

• Schnittdiagramm

Erklärende Variablen sind Daten, die zwei unabhängigen Normalverteilungen folgen.

• Erklärende Variablen $ x ^ {(1)} _ t, x ^ {(2)} _ t $

• Ergebnis $ y $

Infolgedessen scheint der Graph von $ y $ keine Struktur mit nur Zufallszahlen für das menschliche Auge zu haben. Ich möchte bestätigen, ob die latenten Variablen aus diesen Daten gut geschätzt werden können.

3. Modellschätzung mit dlm

Nun, hier ist die Produktion. Ich möchte dlm von R von Python über rpy2 aufrufen, um das zeitvariable Koeffizientenregressionsmodell abzuschätzen.

3-1. Rpy2-Import

Importiere 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()

Importieren Sie die R-Bibliothek, dlm.

dlm = importr("dlm")

3-2. Modellspezifikation des Zustandsraummodells

Hier geben Sie den Typ des Zustandsraummodells an. Da es sich diesmal um ein Regressionsmodell mit erklärenden Variablen handelt, verwenden wir "dlmModReg". "X" ist eine erklärende Variable, und dieses Mal gibt es zwei, "x1" und "x2". Setzen Sie also diese Daten. dV ist die Varianz des Beobachtungsmodellfehlers und repräsentiert die Varianz von $ v_t $ $ \ sigma_v $. dW ist die Verteilung von Systemmodellfehlern. Dieses Mal gibt es drei Systemmodelle, die sich auf $ e ^ {(1)} _t, e ^ {(2)} _t, w_t $ beziehen.

Es entspricht $ \ sigma_ {e} ^ {(1)}, \ sigma_ {e} ^ {(2)} bzw. \ 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])]
  )

Nachdem der Modelltyp festgelegt wurde, stellen Sie den beobachteten Wert "y" ein und ermitteln Sie die Parameter mit der wahrscheinlichsten Methode. Da die Schätzung abhängig vom Anfangswert falsch sein kann, wird sie zweimal unter Verwendung von "dlmMLE" geschätzt. Dies bedeutet, dass das Ergebnis des ersten "dlmMLE" als Anfangswert verwendet wird und das zweite "dlmMLE" ausgeführt wird, um die Stabilität zu verbessern.

3-3. Schätzen der Varianz des Fehlers nach der wahrscheinlichsten Methode

python

#Führen Sie die wahrscheinlichste Schätzung in zwei Schritten durch
%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")

Lassen Sie uns das Ergebnis anzeigen.

python

show_result(fitDlmReg)

Das folgende "par" entspricht den Schätzergebnissen der vier Standardabweichungen von $ \ sigma_v, e ^ {(1)} _t, e ^ {(2)} _t, w_t $, ist jedoch auf "dlmModReg" gesetzt. Manchmal setze ich "np.exp" dazwischen, um zu verhindern, dass der Wert negativ wird. Wenn ich also das Ergebnis "np.exp" und dann "np.sqrt" nehme, ist dies die Standardabweichung.

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

Wenn Sie tatsächlich berechnen, können Sie sehen, dass die beim Erstellen der künstlichen Daten angegebenen Parameter nahezu reproduziert werden.

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]})

3-4 Durchführen von Filtern und Glätten

Suchen Sie nun anhand dieses Ergebnisses den vom Kalman-Filter gefilterten Wert und anschließend den geglätteten Wert aus diesem Filterwert.

python

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

#Glätten
smoothDlmReg = dlm.dlmSmooth(filterDlmReg)

Die erhaltenen Filter- und Glättungsergebnisse werden in die Python-Welt zurückgegeben und als "ndarray" gespeichert.

python

#Erhalten Sie gefilterte Werte
filtered_a = np.array(get_method(filterDlmReg, "a"))

#Erhalten Sie einen geglätteten Wert
smoothed_a = np.array(get_method(smoothDlmReg, "s"))

3-5 Visualisierung der Schätzergebnisse

Unten sehen Sie eine grafische Darstellung der beiden Ergebnisse des Abschnitts "Achsenabschnitt" und des Regressionskoeffizienten $ a ^ {(1)} _ t, a ^ {(2)} _ t $. Der ursprüngliche Wert der künstlichen Daten wird in blau angezeigt, der gefilterte Wert wird in grün angezeigt und der geglättete Wert wird in rot angezeigt, aber ich denke, dass die Bewegung der ursprünglichen Daten erheblich reproduziert werden kann. Natürlich erhält dlm nur drei Daten, "x1, x2, y" und den Typ des Zustandsraummodells, aber der Schnitt und der Regressionskoeffizient können ordnungsgemäß wiederhergestellt werden! (Der Wert ist nur in der Nähe des Anfangswertes weit verbreitet.)

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")

Vergleichen wir abschließend den Wert ohne den Beobachtungsfehler des y-Werts mit dem Wert, der mit y unter Verwendung der geglätteten latenten Variablen berechnet wurde. Sie können auch sehen, dass der wahre y-Wert durch Entfernen des Beobachtungsfehlers erheblich reproduziert werden kann!

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.")

Referenz

[1] Logik des blauen "Zeitvariablen Koeffizientenmodells nach 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] Dieser Code wird hochgeladen (github) https://github.com/matsuken92/Qiita_Contents/blob/master/General/dlm_time_varying_regression.ipynb

[2] Version dieses Codes, die von n Variablen unterstützt werden kann (github) https://github.com/matsuken92/Qiita_Contents/blob/master/General/dlm_time_varying_regression_n-var.ipynb