Benutzerdefiniertes Zustandsraummodell in Python

Einführung

Unter den Zustandsraummodellen gab es nur wenige Websites, auf denen benutzerdefinierte Modelle auf Japanisch erklärt wurden (ich konnte sie nicht finden). Daher werde ich die Ergebnisse meiner eigenen Forschung veröffentlichen.

Für Statistikmodelle (insbesondere benutzerdefinierte Zustandsraummodelle) war es schwierig, das Modell zu definieren, daher hatte ich Schwierigkeiten damit.

Ich würde mich freuen, wenn Sie auf falsche Beschreibungen hinweisen könnten.

Was ist ein Zustandsraummodell?

Beobachtungsgleichung und Zustandsgleichung

Die Struktur des Modells besteht aus einer Beobachtungsgleichung (es gibt auch ein Dokument, das es als Beobachtungsmodell beschreibt) und einer Zustandsgleichung (es gibt auch ein Dokument, das es als Systemmodell beschreibt).

Beobachtungsgleichung: $ y_t = Z_t (x_ {t}) + d_t + \ epsilon_t $ Zustandsgleichung: $ x_ {t + 1} = T_t (x_ {t}) + c_t + R_t \ eta_ {t} $

Vergangenheit, Gegenwart und Zukunft

Obwohl eine detaillierte Erklärung weggelassen wird, nimmt das Zustandsraummodell die folgenden Schätzungen für die Vergangenheit, Gegenwart und Zukunft vor.

Bedeutung des benutzerdefinierten Modells

Unter den Zustandsraummodellen ist die Bedeutung des benutzerdefinierten Modells wie folgt.

Modell verwendet diesmal

Beobachtungsgleichung: 　y_t = x_{t} + \epsilon_t, \hspace{10pt}\epsilon_t\sim N(0,H_t) Zustandsgleichung: 　x_{t} = T_t x_{t-1} + A\cdot\exp\left(\frac{B}{z}\right) + \eta_{t}, \hspace{10pt}\eta_{t}\sim N(0,Q_t)

Der Unterschied zur obigen Formel besteht darin, dass $ Z_t $ in der Beobachtungsgleichung festgelegt ist, $ d_t $ im konstanten Term weggelassen wird und $ c_t $ im konstanten Term in der Zustandsgleichung $ A \ cdot \ ist. Ersetzt durch exp \ left (\ frac {B} {z} \ right) $.

Bibliotheksimport etc.

Hier ist der Python-Code.

import numpy as np
import pandas as pd
import datetime
import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("figure", figsize=(16,8))
plt.rc("font", size=15)

Datenerstellung zum Testen

Erstellen Sie Daten zum Testen. (Wenn Sie es tatsächlich verwenden, ist eine solche Arbeit nicht erforderlich, sodass eine ausführliche Erklärung entfällt.)

def gen_data_for_model1():
    nobs = 1000

    rs = np.random.RandomState(seed=93572)

    Ht = 5
    Tt = 1.001
    A = 0.1
    B = -0.007
    Qt = 0.01
    var_z = 0.1

    et = rs.normal(scale=Ht**0.5, size=nobs)
    z = np.cumsum(rs.normal(size=nobs, scale=var_z**0.05))
    Et = rs.normal(scale=Qt**0.5, size=nobs)

    xt_1 = 50
    x = []
    
    for i in range(nobs):
        xt = Tt * xt_1
        x.append(xt)
        xt_1 = xt
    xt = np.array(x)
    xt = xt + A * np.exp(B/z) + Et
    yt = xt + et
    
    return yt, xt, z

yt, xt, z = gen_data_for_model1()
_ = plt.plot(yt,color = "r")
_ = plt.plot(xt, color="b")

data.png

Speichern Sie Daten in DataFrame. Setzen Sie den Index auf datetime, um Vorhersagen zu treffen

df = pd.DataFrame()
df['y'] = yt
df['x'] = xt
df['z'] = z
st = datetime.datetime.strptime("2001/1/1 0:00", '%Y/%m/%d %H:%M')
date = []
for i in range(1000):
    if i == 0:
        d = st
    dt = d.strftime('%Y/%m/%d %H:%M')
    date.append(dt)
    d += datetime.timedelta(days=1)
df['date'] = date
df['date'] = pd.to_datetime(df['date'] )
df = df.set_index("date")
df

df.png

Modelldefinition

class custom(sm.tsa.statespace.MLEModel):
    param_names = ['T', 'A', 'B', 'Ht', 'Qt'] #Geben Sie die zu schätzenden Parameter an
    start_params = [1., 1., 0., 1., 1] #Anfangswert des zu schätzenden Parameters

    def __init__(self, endog, exog):
        exog = np.squeeze(exog)
        super().__init__(endog, exog=exog, k_states=1, initialization='diffuse')
        self.k_exog = 1
        
        # Z = I,Der Teil, der Zt in der obigen Gleichung entspricht, wird diesmal zur Einheitsmatrix.
        self['design', 0, 0] = 1.
        # R = I,Der Teil, der Rt in der obigen Gleichung entspricht, wird diesmal zur Einheitsmatrix.
        self['selection', 0, 0] = 1.
        
        # c_Stellen Sie t auf Zeit ein
        self['state_intercept'] = np.zeros((1, self.nobs))
        
    def clone(self, endog, exog, **kwargs):
        return self._clone_from_init_kwds(endog, exog=exog, **kwargs)
        
    def transform_params(self, params):
        #Die Varianz muss eine positive Zahl sein, also quadrieren Sie sie einmal
        params[3:] = params[3:]**2
        return params
        
    def untransform_params(self, params):
        #Quadrat 1/Quadrat, um zur ursprünglichen Größe zurückzukehren (Quadratwurzel)
        params[3:] = params[3:]**0.5
        return params
    
    def update(self, params, **kwargs):
        #Es ist eine zu aktualisierende Spezifikation
        params = super().update(params, **kwargs)
        # T = T
        self['transition', 0, 0] = params[0] 
        # c_t = A * exp(B / z_t)
        self['state_intercept', 0, :] = params[1] * np.exp(params[2] / self.exog)       
        # Ht
        self['obs_cov', 0, 0] = params[3]
        # Qt
        self['state_cov', 0, 0] = params[4]

Teilen Sie in Schätzdaten und Vorhersagedaten

df_test = df.iloc[:800,:]
df_train = df.iloc[800:,:]

Lernen, zusammenfassende Bestätigung

mod = custom(df_test['y'], df_test['z'])
res = mod.fit()
print(res.summary())

サマリー.png

Der Coef-Teil ist der geschätzte Wert, liegt jedoch nahe am numerischen Wert der zu Beginn erstellten Daten.

Speichern Sie das Schätzergebnis in DataFrame

ss = pd.DataFrame(res.smoothed_state.T, columns=['x'], index=df_test.index)

Schätzung und Prognose

predict = res.get_prediction()
forecast = res.get_forecast(df.index[-1], exog = df_train['z'].values)

Bei der Vorhersage müssen Sie angeben, bis zu welchem Punkt Sie eine Vorhersage treffen möchten. (Diesmal angegeben durch df.index [-1]) Wenn exogene Variablen im Modell verwendet werden, müssen die Daten zum Zeitpunkt der Vorhersage mit dem Argument "exog =" angegeben werden. Dieses Mal sind alles Dummy-Daten, aber wenn Sie sie tatsächlich für die Vorhersage verwenden möchten, müssen Sie Dummy-Daten aus früheren Daten erstellen. Daten mit ungleichmäßigem Datenabstand können nicht zur Vorhersage verwendet werden. (Gleiches gilt für SARIMAX, das auch exogene Variablen verwendet.)

Visualisierung

fig, ax = plt.subplots()

y = df['y']

y.plot(ax=ax,label='y')
predict.predicted_mean.plot(label='x')
predict_ci = predict.conf_int(alpha=0.05)
predict_index = predict_ci.index
ax.fill_between(predict_index[2:], predict_ci.iloc[2:, 0], predict_ci.iloc[2:, 1], alpha=0.1)

forecast.predicted_mean.plot(ax=ax, style='r', label='forecast')
forecast_ci = forecast.conf_int()
forecast_index = forecast_ci.index
ax.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], alpha=0.1)

legend = ax.legend(loc='best');
fig.savefig('custom_statespace.png')

可視化.png

Der hellblaue Teil ist das geschätzte 95% -Konfidenzintervall und der rosa Teil ist das vorhergesagte 95% -Konfidenzintervall.

Verweise

Dieser Artikel wurde unter Bezugnahme auf die folgenden Informationen verfasst.

Vielen Dank an Chad Fulton, den Entwickler von Statistikmodellen!

Empfohlene Bücher für Zustandsraummodelle

Wir werden auch andere empfohlene Bücher als die oben aufgeführten vorstellen. Sie sind in der Reihenfolge ihrer relativen Lesbarkeit aufgeführt.

Bücher, die ich gekauft, aber noch nicht gelesen habe. .. ..

Versionsgeschichte

  1. Oktober 2020: Erstausgabe

Recommended Posts

Benutzerdefiniertes Zustandsraummodell in Python
Benutzerdefinierte Sortierung in Python3
Visualisieren Sie das Keras-Modell mit Python 3.5
Allgemeines Gaußsches Zustandsraummodell in Python
Lernen Sie das Entwurfsmuster "State" in Python
Implementieren Sie ein benutzerdefiniertes Benutzermodell in Django
Python in der Optimierung
CURL in Python
Metaprogrammierung mit Python
Python 3.3 mit Anaconda
Geokodierung in Python
Metaanalyse in Python
Unittest in Python
Epoche in Python
Zwietracht in Python
Deutsch in Python
DCI in Python
Quicksort in Python
nCr in Python
N-Gramm in Python
Programmieren mit Python
Plink in Python
Konstante in Python
FizzBuzz in Python
Schritt AIC in Python
CSV in Python
Reverse Assembler mit Python
Reflexion in Python
Implementierung des Partikelfilters durch Python und Anwendung auf das Zustandsraummodell
Konstante in Python
nCr in Python.
Format in Python
Scons in Python 3
Puyopuyo in Python
Python in Virtualenv
PPAP in Python
Quad-Tree in Python
Chemie mit Python
Hashbar in Python
DirectLiNGAM in Python
LiNGAM in Python
In Python reduzieren
In Python flach drücken
Ich habe versucht, TOPIC MODEL in Python zu implementieren
Erstellen Sie in Python ein einfaches Momentum-Investmentmodell
Verwenden Sie eine benutzerdefinierte Fehlerseite mit Python / Tornado
Übertragungsfunktion / Zustandsraummodell der RLC-Serienschaltung und Simulation von Python
Täglicher AtCoder # 36 mit Python
Clustertext in Python
Täglicher AtCoder # 32 in Python
Täglicher AtCoder # 6 in Python
Täglicher AtCoder # 18 in Python
Bearbeiten Sie Schriftarten in Python
Singleton-Muster in Python
Dateioperationen in Python
Lesen Sie DXF mit Python
Täglicher AtCoder # 53 in Python
Übertragungsfunktion / Zustandsraummodell des Feder- / Masse- / Dämpfersystems und Simulation von Python