Modèle d'espace d'états personnalisé en Python

introduction

Parmi les modèles d'espace d'état, il y avait peu de sites qui expliquaient les modèles personnalisés en japonais (je ne pouvais pas les trouver), je publierai donc les résultats de mes propres recherches.

Pour les modèles de statistiques (en particulier les modèles d'espace d'état personnalisés), il était difficile de définir le modèle, donc j'ai eu du mal à le faire.

Je vous serais reconnaissant si vous pouviez signaler des descriptions incorrectes.

Qu'est-ce qu'un modèle d'espace d'états?

Équation d'observation et équation d'état

La structure du modèle se compose d'une équation d'observation (il existe également un document qui la décrit comme un modèle d'observation) et d'une équation d'état (il existe également un document qui la décrit comme un modèle système).

Équation d'observation: $ y_t = Z_t (x_ {t}) + d_t + \ epsilon_t $ Équation d'état: $ x_ {t + 1} = T_t (x_ {t}) + c_t + R_t \ eta_ {t} $

Passé, présent et futur

Bien qu'une explication détaillée soit omise, le modèle d'espace d'états effectue les estimations suivantes pour le passé, le présent et le futur.

Importance du modèle personnalisé

Parmi les modèles d'espace d'états, la signification du modèle personnalisé est la suivante.

Modèle utilisé cette fois

Équation d'observation: 　y_t = x_{t} + \epsilon_t, \hspace{10pt}\epsilon_t\sim N(0,H_t) Équation d'état: 　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)

La différence avec la formule ci-dessus est que $ Z_t $ est fixé dans l'équation d'observation, $ d_t $ dans le terme constant est omis et $ c_t $ dans le terme constant dans l'équation d'état est $ A \ cdot . Remplacé par exp \ left (\ frac {B} {z} \ right) $.

Importation de bibliothèque, etc.

Voici maintenant le code Python.

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)

Création de données pour les tests

Créez des données pour les tests. (Lors de son utilisation réelle, un tel travail est inutile, une explication détaillée est donc omise.)

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

Stockez les données dans DataFrame. Définir l'index sur datetime pour faire des prédictions

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

Définition du modèle

class custom(sm.tsa.statespace.MLEModel):
    param_names = ['T', 'A', 'B', 'Ht', 'Qt'] #Spécifiez les paramètres à estimer
    start_params = [1., 1., 0., 1., 1] #Valeur initiale du paramètre à estimer

    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,La partie correspondant à Zt dans l'équation ci-dessus devient cette fois la matrice unitaire.
        self['design', 0, 0] = 1.
        # R = I,La partie correspondant à Rt dans l'équation ci-dessus devient cette fois la matrice unitaire.
        self['selection', 0, 0] = 1.
        
        # c_Régler t à l'heure
        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):
        #La variance doit être un nombre positif, alors mettez-la au carré une fois
        params[3:] = params[3:]**2
        return params
        
    def untransform_params(self, params):
        #Au carré 1/Carré pour revenir à la taille d'origine (racine carrée)
        params[3:] = params[3:]**0.5
        return params
    
    def update(self, params, **kwargs):
        #C'est une spécification à mettre à jour
        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]

Diviser en données d'estimation et données de prédiction

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

Apprentissage, confirmation sommaire

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

サマリー.png

La partie coef est la valeur estimée, mais elle est proche de la valeur numérique des données créées au début.

Stocker le résultat de l'estimation dans DataFrame

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

Estimation et prévision

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

Lors de la prédiction, vous devez spécifier à quel point vous souhaitez prédire. (Cette heure est spécifiée par df.index [-1]) De plus, lors de l'utilisation de variables exogènes dans le modèle, il est nécessaire de spécifier les données avec l'argument "exog =" au moment de la prédiction. Cette fois, tout est des données factices, mais si vous souhaitez réellement les utiliser dans la prédiction, vous devez créer des données factices à partir de données passées. Les données dont l'espacement des données est irrégulier ne peuvent pas être utilisées pour la prédiction. (Il en va de même pour SARIMAX, qui utilise également des variables exogènes.)

Visualisation

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

La partie bleu clair est l'intervalle de confiance estimé à 95% et la partie rose est l'intervalle de confiance prévu à 95%.

Les références

Cet article a été rédigé en référence aux informations suivantes.

Merci à Chad Fulton, le développeur de statsmodels!

Livres recommandés pour les modèles d'espace d'état

Nous présenterons également des livres recommandés autres que ceux énumérés ci-dessus. Ils sont répertoriés par ordre de lisibilité relative.

Livres que j'ai achetés mais que je n'ai pas encore lus. .. ..

Historique des révisions

3 octobre 2020: Première édition

Recommended Posts

Modèle d'espace d'états personnalisé en Python
Tri personnalisé en Python3
Visualiser le modèle Keras avec Python 3.5
Modèle d'espace d'état gaussien général en Python
Apprenez le modèle de conception "État" en Python
Implémenter un modèle utilisateur personnalisé dans Django
Python en optimisation
CURL en Python
Métaprogrammation avec Python
Python 3.3 avec Anaconda
Géocodage en python
Méta-analyse en Python
Unittest en Python
Époque en Python
Discord en Python
Allemand en Python
DCI en Python
tri rapide en python
nCr en python
N-Gram en Python
Programmation avec Python
Plink en Python
Constante en Python
FizzBuzz en Python
Étape AIC en Python
CSV en Python
Assemblage inversé avec Python
Réflexion en Python
Implémentation du filtre à particules par Python et application au modèle d'espace d'états
Constante en Python
nCr en Python.
format en python
Scons en Python 3
Puyopuyo en python
python dans virtualenv
PPAP en Python
Quad-tree en Python
Chimie avec Python
Hashable en Python
DirectLiNGAM en Python
LiNGAM en Python
Aplatir en Python
Aplatir en python
J'ai essayé d'implémenter TOPIC MODEL en Python
Créer un modèle d'investissement dynamique simple en Python
Utilisez une page d'erreur personnalisée avec python / tornado
Fonction de transmission / modèle d'espace d'état du circuit série RLC et simulation par Python
AtCoder # 36 quotidien avec Python
Texte de cluster en Python
Daily AtCoder # 32 en Python
Daily AtCoder # 6 en Python
Daily AtCoder # 18 en Python
Modifier les polices en Python
Motif singleton en Python
Opérations sur les fichiers en Python
Lire DXF avec python
Daily AtCoder # 53 en Python
Fonction de transmission / modèle d'espace d'état du système ressort / masse / amortisseur et simulation par Python