Allgemeines Gaußsches Zustandsraummodell in Python

Einführung

Das Zustandsraummodell ist eine Methode zur Zeitreihenanalyse. Es ist sehr gut interpretierbar, da es modelliert werden kann, indem Änderungen im Laufe der Zeit in mehrere Komponenten zerlegt werden. In den letzten Jahren wurden mehrere Bücher veröffentlicht, die das Zustandsraummodell auf leicht verständliche Weise erklären.

In diesem Artikel nicht erklärt

Siehe die Nachschlagewerke und Kommentarartikel unten.

Was in diesem Artikel zu erklären

Referenzen für Zustandsraummodelle

  1. [Grundlagen der Zeitreihenanalyse und des Zustandsraummodells: Theorie und Implementierung mit R und Stan gelernt](https://www.amazon.co.jp/%E6%99%82%E7%B3%BB%E5%] 88% 97% E5% 88% 86% E6% 9E% 90% E3% 81% A8% E7% 8A% B6% E6% 85% 8B% E7% A9% BA% E9% 96% 93% E3% 83% A2% E3% 83% 87% E3% 83% AB% E3% 81% AE% E5% 9F% BA% E7% A4% 8E-R% E3% 81% A8Stan% E3% 81% A7% E5% AD% A6% E3% 81% B6% E7% 90% 86% E8% AB% 96% E3% 81% A8% E5% AE% 9F% E8% A3% 85-% E9% A6% AC% E5% A0% B4 -% E7% 9C% 9F% E5% 93% 89 / dp / 4903814874 / ref = sr_1_1? __ mk_ja_JP =% E3% 82% AB% E3% 82% BF% E3% 82% AB% E3% 83% 8A & Schlüsselwörter =% E7% 8A% B6% E6% 85% 8B% E7% A9% BA% E9% 96% 93% E3% 83% A2% E3% 83% 87% E3% 83% AB & qid = 1554681669 & s = Gateway & sr = 8-1) Von den drei vorgestellten Büchern hat es die wenigsten Formeln und ist im Klartext geschrieben. Empfohlen für diejenigen, die das Bild des Zustandsraummodells vorerst erfassen und in R implementieren möchten. Es wird hauptsächlich von KFAS von R implementiert, berührt aber auch dlm.

  2. [Zeitreihenanalyse aus dem in R praktizierten Basics-Kalman-Filter, MCMC, Partikelfilter](https://www.amazon.co.jp/%E5%9F%BA%E7%A4%8E%E3%81] % 8B% E3% 82% 89% E3% 82% 8F% E3% 81% 8B% E3% 82% 8B% E6% 99% 82% E7% B3% BB% E5% 88% 97% E5% 88% 86 % E6% 9E% 90-% E2% 80% 95R% E3% 81% A7% E5% AE% 9F% E8% B7% B5% E3% 81% 99% E3% 82% 8B% E3% 82% AB% E3% 83% AB% E3% 83% 9E% E3% 83% B3% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF% E3% 83% BBMCMC% E3% 83% BB% E7% B2% 92% E5% AD% 90% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF% E3% 83% BC-Data-Science- Library / dp / 4774196460 / ref = pd_sbs_0_1 / 356-7900534-8904633? _encoding = UTF8 & pd_rd_i = 4774196460 & pd_rd_r = 6cba9d2a-5991-11e9-b4dc-811bcb7ea7c6 & pd_rd_w = hpZLR & pd_rd_wg = 8TNR0 & pf_rd_p = ad2ea29d-ea11-483c-9db2-6b5875bb9b73 & pf_rd_r = PQSTK8PGCF3HMAA4XW95 & psc = 1 & refRID = PQSTK8PGCF3HMAA4XW95) Ein herzliches Buch. Es werden auch Wiener-Filter angesprochen, die in anderen Büchern nicht erwähnt werden. Sie können auch lernen, wie Sie ein angewendetes Modell mit viel Code implementieren. Verwendet die DLM-Bibliothek von R.

  3. [Kalman Filter-Time Series Prediction und State Space Model Using R-](https://www.amazon.co.jp/%E3%82%AB%E3%83%AB%E3%83%9E% E3% 83% B3% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF-% E2% 80% 95R% E3% 82% 92% E4% BD% BF% E3 % 81% A3% E3% 81% 9F% E6% 99% 82% E7% B3% BB% E5% 88% 97% E4% BA% 88% E6% B8% AC% E3% 81% A8% E7% 8A % B6% E6% 85% 8B% E7% A9% BA% E9% 96% 93% E3% 83% A2% E3% 83% 87% E3% 83% AB% E2% 80% 95-% E7% B5% B1% E8% A8% 88% E5% AD% A6One-Point-2 / dp / 4320112539 / ref = pd_sim_0_3 / 356-7900534-8904633? rhkfr & pf_rd_p = b88353e4-7ed3-4da1-bc65-341dfa3a88ce & pf_rd_r = 7Q1J7M8DSYWV6M0ZDF7G & psc = 1 & refRID = 7Q1J7M8DSYWV6M0ZDF7G) Das höflichste Buch mit Erklärungen basierend auf mathematischen Formeln. Wie der Titel schon sagt, wird es am meisten empfohlen, um den Kalman-Filter zu verstehen. Es gibt auch eine Erklärung des Partikelfilters. Der Beispielcode von KFAS ist ebenfalls enthalten.

Alle enthielten Erklärungen mit mathematischen Formeln und R-Code, was sehr hilfreich war. Andererseits kann in letzter Zeit eine ** Implementierung in Python ** erforderlich sein.

Eine Bibliothek mit dem Namen ** statsmodels ** wird in Python für Zustandsraummodelle verwendet, aber es ist schwierig zu verstehen, wie sie verwendet wird, da die Erklärung der offiziellen Referenz zu klein ist. Die Grundlagen von Zustandsraummodellen, die Statistikmodelle verwenden, finden Sie in diesem Artikel.

Dieser Artikel beschreibt eine etwas fortgeschrittenere Verwendung von Statistikmodellen.

Allgemeines Gaußsches Zustandsraummodell nach Statistikmodellen

Das allgemeine Gaußsche Zustandsraummodell wird durch die folgende Gleichung ausgedrückt.

Das Modell auf lokaler Ebene, das das grundlegendste Zustandsraummodell ist, ist ein Modell, bei dem alle Koeffizientenmatrizen Einheitsmatrizen sind. Durch geeignete Angabe der Koeffizientenmatrix können verschiedene Modelle wie ein Trendkomponentenmodell und ein periodisches Komponentenmodell ausgedrückt werden. Mit KFAS in R und Statistikmodellen in Python können Sie Modelle erstellen, indem Sie die Koeffizientenmatrix direkt angeben.

Es gibt jedoch eine einfachere Möglichkeit, das Grundmodell zu implementieren.

UnobservedComponents-Klasse

Bei der Implementierung eines Modells, das Trends, periodische Komponenten, Regressionskomponenten usw. in Statistikmodellen enthält, im Grunde genommen ** Nicht beobachtete Komponenten ** Verwenden Sie eine Klasse mit dem Namen .structural.UnobservedComponents.html # statsmodels.tsa.statespace.structural.UnobservedComponents).

Bei UnobservedComponents ist der Freiheitsgrad bei der Modellierung jedoch gering, z. B. die Unfähigkeit, ein Regressionskomponentenmodell unter Verwendung zeitvariabler Koeffizienten festzulegen, die Unfähigkeit, mehrere periodische Komponenten einzubeziehen, und die Unfähigkeit, den Zustand mit bekannten Werten zu initialisieren.

Der Vorteil des Kalman-Filters besteht beispielsweise darin, dass er ** sequentielle Vorhersagen ** treffen kann und in der Praxis häufig ** Online-Vorhersagen ** treffen sollte. In diesem Fall wird der Kalman-Filter zuerst für einen bestimmten Zeitraum für die Daten ausgeführt, und dann werden die neu erhaltenen Daten erneut gefiltert. Zu diesem Zeitpunkt werden der Anfangswert des Zustands, der für die nächste Filterung verwendet werden soll, und die Vorhersagefehlerdispersion basierend auf dem Zustandsvorhersagewert zum Endzeitpunkt des vorherigen Kalman-Filters bestimmt.

Unobserved Components nähert sich jedoch standardmäßig der abgelenkten Initialisierung an (einfach weil man denkt, dass es keine Hinweise auf den Anfangszustand gibt) und ** kann nicht mit bekannten Werten initialisiert werden **. In diesen Fällen müssen Sie das MLE-Modell verwenden.

MLEModel-Klasse

Wenn Sie Ihr eigenes Modell festlegen möchten, lesen Sie bitte ** MLE-Modell **, wie unter Offizielle Referenz erläutert. https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.mlemodel.MLEModel.html?highlight=mlemodels)というクラスを使えば良い。 Auf der offiziellen Website finden Sie Erläuterung zur Implementierung des Trendmodells des lokalen Leitungssystems. Dieses Modell kann jedoch mit nicht beobachteten Komponenten implementiert werden.

Hier werden wir die Implementierung eines Modells vorstellen, das die folgenden zwei Punkte erfüllt, die von nicht beobachteten Komponenten nicht realisiert werden können.

  1. Regressionskomponentenmodell mit zeitlich variierendem Regressionskoeffizienten
  2. Die Initialisierungsmethode kann ausgewählt werden

CustomizeMLEmodel.py


import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

"""
Regression Model
"""
class StateSpaceRegression(sm.tsa.statespace.MLEModel):
    def __init__(self, endog, exog, initialization, initial_state=None, initial_state_cov=None):
        # Model order
        # k_state:Zustandsdimension
        # k_posdef:Prozessfehler Dispersionsmatrixabmessungen
        k_states = k_posdef = 2

        # Initialize the statespace
        super(StateSpaceRegression, self).__init__(
            endog, exog=exog, k_states=k_states, k_posdef=k_posdef,
            initialization=initialization, initial_state=initial_state, initial_state_cov=initial_state_cov,
            loglikelihood_burn=k_states
        )

        # Initialize the matrices
        #Die erklärenden Variablen werden durch eine zeitlich variierende Entwurfsmatrix dargestellt.
        n_obs = len(exog)
        self.ssm['design'] = np.vstack((np.repeat(1, n_obs), exog))[np.newaxis, :, :]
        
        self.ssm['transition'] = self.ssm['selection'] = np.eye(k_states)

        #Selbst mit Update.Zustand zu ssm_Wird verwendet, um den Wert von cov festzulegen
        self._state_cov_idx = ('state_cov',) + np.diag_indices(k_posdef)

    @property
    def param_names(self):
        return ['sigma2.measurement', 'sigma2.intercept', 'sigma2.coefficient']

    @property
    #Anfangswert bei der Schätzung des wahrscheinlichsten Parameters
    def start_params(self):
        return [np.std(self.endog)]*3

    def transform_params(self, unconstrained):
        return unconstrained**2

    def untransform_params(self, constrained):
        return constrained**0.5

    def update(self, params, *args, **kwargs):
        params = super(StateSpaceRegression, self).update(params, *args, **kwargs)
        
        # Observation covariance
        self.ssm['obs_cov',0,0] = params[0]

        # State covariance
        self.ssm[self._state_cov_idx] = params[1:]

Die Punkte sind wie folgt.

  1. k_state, k_posdef: Die Dimensionen des Zustandsvektors und der Prozessfehlerdispersionsmatrix. Da die Dimension der Zielvariablen (beobachteter Wert) aus den Daten abgeleitet wird, ist es nicht erforderlich, sie anzugeben.
  2. Update: Wird verwendet, wenn der wahrscheinlichste Parameter mit der Anpassungsmethode geschätzt wird. Nimmt Parameter und legt die entsprechende Zustandsraummatrix fest.
  3. Statespace-Matrizen: Standardmäßig sind alle Matrizen Nullmatrizen. Beachten Sie, dass es sich nicht um eine Einheitsmatrix handelt. Wenn es sich mit der Zeit ändert, muss es mit der Aktualisierungsmethode festgelegt werden. Es ist in self.ssm gespeichert.
  4. Startparameter: Anfangswerte für die wahrscheinlichste Schätzung von Verteilungsparametern
  5. Initialisierung: Sie müssen den Mittelwert und die Varianz der Anfangsverteilung des Zustandsvektors festlegen. Wenn die anfängliche Verteilung bekannt ist, wird initialize_known verwendet, und wenn sie unbekannt ist, wird initialize_approximate_diffuse (ungefähre diffuse Initialisierung) verwendet.
  6. transformieren: Standardmäßig wird die Parameterschätzung ohne Einschränkungen durchgeführt. Durch Festlegen von transform_params kann sie jedoch in eingeschränkte Parameter wie Verteilungen transformiert werden, die positive Werte annehmen.
  7. param_names: Sie können die Parameter benennen, auf die Sie schließen möchten.

Anwendung auf reale Daten

Hier sind die Ergebnisse von Experimenten mit Beispieldaten, die auf der Support-Website in der dritten Referenz oben verfügbar sind.

Das erste Diagramm ist der gemessene Wert. $ Y $ (blau) wird durch ein Regressionskomponentenmodell des zeitvariablen Koeffizienten mit $ x $ (orange) als erklärende Variable vorhergesagt. Das zweite Diagramm repräsentiert die Pegelkomponente und das dritte Diagramm repräsentiert den Regressionskoeffizienten. Im Fall dieser Daten ist ersichtlich, dass der Regressionskoeffizient tendenziell gegen die zweite Hälfte des Zeitraums ansteigt.

nintendo.png

Klicken Sie hier für den Code.

nintendo = pd.read_csv('data/NINTENDO.csv', index_col=0, parse_dates=[0])
nikkei = pd.read_csv('data/NIKKEI225.csv', index_col=0, parse_dates=[0])

def scaler(s):
    return (s - s.mean()) / s.std()

endog, exog = scaler(nintendo.Close), scaler(nikkei.Close)

#Modell des zeitvariablen Regressionskoeffizienten
mod_reg = StateSpaceRegression(endog, exog, initialization='approximate_diffuse')

res_mod_reg = mod_reg.fit()
res_mod_reg.summary()

#Ein-Jahres-Prognose des Staates
mu_pred, beta_pred  = res_mod_reg.predicted_state
#Geglätteter Zustand
mu_smoothed, beta_smoothed = res_mod_reg.smoothed_state


fig, axes = plt.subplots(3, 1, figsize=(8, 8), sharex=True)
endog.plot(ax=axes[0], label='y')
exog.plot(ax=axes[0], label='x')
axes[0].set_title('observed variables')
axes[0].legend()

axes[1].plot(endog.index, mu_smoothed)
axes[1].set_title('level component')

axes[2].plot(endog.index, beta_smoothed)
axes[2].set_title('regression coefficient')

Bonus: TensorFlow-Wahrscheinlichkeitsoption

Bei Verwendung von Statistikmodellen werden die Parameter mit der wahrscheinlichsten Methode + Kalman-Filter punktgeschätzt. Wenn jedoch eine basianische Modellierung möglich ist, kann ein flexibles Modell erstellt und geschätzt werden, und die Sicherheit der Vorhersage kann auf natürliche Weise quantifiziert werden. Es gibt. Es scheint Standard zu sein, Stan zur Implementierung des Bayes'schen Zustandsraummodells zu verwenden, aber es erfordert zwangsläufig Codierungsaufwand. Auch wenn es Pystan gibt, habe ich das Gefühl, dass es oft in R verwendet wird.

Bei der Implementierung des Bayes'schen Zustandsraummodells in Python ist das TensorFlow Probability (TFP) -Sts-Modul eine neue Option. Siehe Offizieller Artikel.

Bis Dezember 2019 befindet sich TFP jedoch noch in der Entwicklung, sodass es noch nicht stabil ist. Erwartungen für die Zukunft!

Recommended Posts

Allgemeines Gaußsches Zustandsraummodell in Python
Benutzerdefiniertes Zustandsraummodell in Python
Allgemeine Relativitätstheorie in Python: Einführung
Zeitreihenanalyse durch allgemeines Gaußsches Zustandsraummodell unter Verwendung von Python [Implementierungsbeispiel unter Berücksichtigung von Extrinsik und Saisonalität]
Quadtree in Python --2
Python in der Optimierung
CURL in Python
Metaprogrammierung mit Python
Python 3.3 mit Anaconda
Geokodierung in Python
SendKeys in Python
Metaanalyse in Python
Unittest in Python
Epoche 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
SQLite in Python
Schritt AIC in Python
LINE-Bot [0] in Python
CSV in Python
Reverse Assembler mit Python
Reflexion in Python
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
Reflexion in Python
Chemie mit Python
Hashbar in Python
DirectLiNGAM in Python
LiNGAM in Python
In Python reduzieren
In Python flach drücken
[Python] Clustering mit einem unendlich gemischten Gaußschen Modell
Ich habe versucht, TOPIC MODEL in Python zu implementieren
Erstellen Sie in Python ein einfaches Momentum-Investmentmodell
Sortierte Liste in Python
Täglicher AtCoder # 36 mit Python
Clustertext in Python
AtCoder # 2 jeden Tag mit Python
[Python] Implementierung von Clustering mit einem gemischten Gaußschen Modell
Täglicher AtCoder # 32 in Python
Täglicher AtCoder # 6 in Python
Täglicher AtCoder # 18 in Python
Bearbeiten Sie Schriftarten in Python
Dateioperationen in Python
Lesen Sie DXF mit Python