Markov-Schaltmodell von Python

1. Wirtschaftsdaten "Regime"

Wenn eine plötzliche und große Veränderung in der Wirtschaft eintritt, wie der Zusammenbruch der IT-Blase oder der Lehman-Schock, kann davon ausgegangen werden, dass sich ein Staat im Hintergrund grundlegend verändert hat. Ein solches grundlegendes Zustandsänderungsmodell wird als "Regime Switching Model" bezeichnet. Das Modellieren eines Regimewechsels kann schwierig sein, wenn der geänderte Zustand nicht beobachtbar ist. Die Modellierung ist jedoch möglich, indem die Wahrscheinlichkeitsverteilung des Regimes mit Parametern angegeben und der Übergang zwischen Regimen als Markov-Kette behandelt wird. Die Parameter, die die Wahrscheinlichkeit maximieren (Mittelwert / Standardabweichung der Normalverteilung, Übergangswahrscheinlichkeit des Regimes), können von MCMC erhalten werden.

2. Erläuterung des Markov-Schaltmodells

Lineare Regression mit Regime

Betrachten Sie das Regime-Switching-Modell des folgenden linearen Regressionsmodells.

\begin{eqnarray}
y\left( t \right) &=& x\left( t \right)\beta + \epsilon\left( t \right)\\
\epsilon\left( t \right)&\sim& N\left(0,\sigma^2\right)
\end{eqnarray}

Die $ 1 $ bis $ n $ Regime zu jedem Zeitpunkt $ t $ werden wie folgt ausgedrückt.

\begin{eqnarray}
I\left(t\right) \in { 1,2,\cdots ,n}
\end{eqnarray}

Wenn das Regime beobachtet werden kann, wird die Regressionsgleichung wie folgt ausgedrückt, wobei das Regime als Index verwendet wird.

\begin{eqnarray}
y\left( t \right) &=& x\left( t \right)\beta_{I\left( t\right)} + \epsilon\left( t \right)\\
\epsilon\left( t \right)&\sim& N\left(0,\sigma^2_{I\left( t\right)}\right)
\end{eqnarray}

Zu diesem Zeitpunkt ist $ y \ left (t \ right) $ der Durchschnitt von $ x \ left (t \ right) \ beta_ {I \ left (t \ right)} $, verteilt auf $ \ sigma ^ 2_ {I \ left ( Es folgt eine Normalverteilung von t \ right)} $. Wenden Sie die wahrscheinlichste Methode an, die die logarithmische Wahrscheinlichkeit maximiert. Das heißt, der Parameter, der die folgende Wahrscheinlichkeitsfunktion maximiert, wird erhalten.

\begin{eqnarray}
\ln{\mathcal{L}}&=&\sum_{t=1}^{T}\ln{\{f\left( y\left( t\right)|I\left( t\right)\right)\} }\\
f\left( y\left( t\right)|I\left( t\right)\right)&=&\frac{1}{\sqrt{2\pi\sigma^2_{I\left( t\right)}}}\exp\left( -\frac{\{y\left( t\right)-x\left( t\right)\beta_{I\left( t \right)}\}^2}{2\sigma^2_{I\left( t\right)}}\right)
\end{eqnarray}

Wenn das Regime nicht direkt beobachtet werden kann, wird auch angenommen, dass $ I \ left (t \ right) $ einer bedingten Wahrscheinlichkeitsverteilung in der Vergangenheit folgt. Sei $ \ mathcal {F} \ left (t \ right) $ die Information, die zum Zeitpunkt $ t $ erhalten wird. Beachten Sie, dass $ f \ left (t \ right) $ nicht von $ \ mathcal {F} \ left (t \ right) $ abhängig ist.

\begin{eqnarray}
f\left( y\left( t\right) | \mathcal{F\left( t-1\right)}\right) &=& \sum_{I\left( t\right)=1}^{n}f\left( y\left( t\right),I\left( t\right) |\mathcal{F}\left( t-1\right)\right)\\
&=& \sum_{I\left( t\right)=1}^{n}f\left( y\left( t\right)|I\left( t\right),\mathcal{F}\left( t-1\right)\right)f\left( I\left( t\right)|\mathcal{F}\left( t-1\right)\right)\\
&=& \sum_{I\left( t\right)=1}^{n}\frac{1}{\sqrt{2\pi\sigma^2_{I\left( t\right)}}}\exp\left( -\frac{\{y\left( t\right)-x\left( t\right)\beta_{I\left( t \right)}\}^2}{2\sigma^2_{I\left( t\right)}}\right)
f\left( I\left( t\right)|\mathcal{F}\left( t-1\right)\right)
\end{eqnarray}

Daher ist die zu maximierende Protokollwahrscheinlichkeit wie folgt.

\begin{eqnarray}
\ln{\mathcal{L}}&=&\sum_{t=1}^{T}\ln{\left\{\sum_{I\left( t\right)=1}^{n}f\left( y\left( t\right)|I\left( t\right),\mathcal{F}\left( t-1\right)\right)f\left( I\left( t\right)|\mathcal{F}\left( t-1\right)\right)\right\} }\\

f\left( y\left( t\right)|I\left( t\right),\mathcal{F}\left( t-1\right)\right)&=&\frac{1}{\sqrt{2\pi\sigma^2_{I\left( t\right)}}}\exp\left( -\frac{\{y\left( t\right)-x\left( t\right)\beta_{I\left( t \right)}\}^2}{2\sigma^2_{I\left( t\right)}}\right)
\end{eqnarray}

Markov-Kette

Wahrscheinlichkeitsverteilung von $ y \ left (t \ right) $ durch Angabe von $ f \ left (I \ left (t \ right) | \ mathcal {F} \ left (t-1 \ right) \ right) $ Ist entschieden. Angenommen, $ I \ left (t \ right) $ erfüllt die Markov-Eigenschaft erster Ordnung.

\begin{eqnarray}
P\left( I\left( t\right)=i|\mathcal{F}\left( t-1\right)\right)
&=&\sum_{j=1}^{n}P\left( I\left( t\right)=i,I\left( t-1\right)=j|\mathcal{F}\left( t-1\right)\right)\\
&=&\sum_{j=1}^{n}P\left( I\left( t\right)=i|I\left( t-1\right)=j\right)P\left( I\left( t-1\right)=j|\mathcal{F}\left( t-1\right)\right)
\end{eqnarray}

Übergangswahrscheinlichkeit von Zustand j zu Zustand i $ P \ links (I \ links (t \ rechts) = i | I \ links (t-1 \ rechts) = j \ rechts) $ zu $ p_ {i, j} $ Betrachten Sie die folgende Übergangswahrscheinlichkeitsmatrix.

\begin{eqnarray}
P &=&
\begin{pmatrix}
p_{1,1}&p_{1,2}&\cdots&p_{1,n}\\
p_{2,1}&\ddots&&\vdots\\
\vdots&&\ddots&\vdots\\
p_{n,1}&\cdots&\cdots&p_{n,n}
\end{pmatrix}
\end{eqnarray}

Bedingte Wahrscheinlichkeit des Regimes $ P \ left (I \ left (t \ right) = j | \ mathcal {F} \ left (t \ right) \ right) bei gegebener Information bis zum Zeitpunkt $ t $ $ Wird durch das Bayes'sche Update berechnet.

\begin{eqnarray}
P\left( I\left( t\right)=i|\mathcal{F}\left( t\right)\right)
&=&P\left( I\left( t\right)=i|\mathcal{F}\left( t-1\right), y\left( t\right)\right)\\
&=&\frac{f\left( I\left( t\right) =i,y\left( t\right)|\mathcal{F}\left( t-1\right)\right)}{f\left(y\left( t\right)|\mathcal{F}\left( t-1\right)\right)}\\
&=&\frac{f\left(y\left( t\right)|I\left( t\right)=i,\mathcal{F}\left( t-1\right)\right)f\left(I\left( t\right)=i|\mathcal{F}\left( t-1\right)\right)}{\sum_{i=1}^{n}f\left(y\left( t\right)|I\left( t\right)=i,\mathcal{F}\left( t-1\right)\right)f\left(I\left( t\right)=i|\mathcal{F}\left( t-1\right)\right)}\\
&=&\frac{\sum_{j=1}^{n}f\left(y\left( t\right)|I\left( t\right)=j,\mathcal{F}\left( t-1\right)\right)p_{i,j}f\left(I\left( t-1\right)=j|\mathcal{F}\left( t-1\right)\right)}{\sum_{i=1}^{n}\sum_{j=1}^{n}f\left(y\left( t\right)|I\left( t\right)=j,\mathcal{F}\left( t-1\right)\right)p_{i,j}f\left(I\left( t-1\right)=i|\mathcal{F}\left( t-1\right)\right)}
\end{eqnarray}

Die Anfangswahrscheinlichkeit ist die stationäre Wahrscheinlichkeit der Markov-Kette. Es ist bekannt, dass die Übergangswahrscheinlichkeitsmatrix $ P $, die Einheitsmatrix $ I $ und der Vektor $ 1 $ mit allen Elementen von $ 1 $ verwendet werden und die stationäre Wahrscheinlichkeit $ \ pi ^ {\ ast} $ wie folgt erhalten werden kann.

\begin{eqnarray}
A &=& 
\begin{pmatrix}
I-P\\
1^{\mathrm{T}}
\end{pmatrix}\\

\pi^{\ast}&=&\left(A^{\mathrm{T}}A\right)^{-1}A^{\mathrm{T}}M.+1. Reihe
\end{eqnarray}

Markov-Ketten-Monte-Carlo-Methode

Die Markov-Ketten-Monte-Carlo-Methode (MCMC) wird als Methode zum Erhalten der Parameter verwendet, die die logarithmische Wahrscheinlichkeit maximieren. Die Erklärung wird hier weggelassen.

Die vorgeschlagenen Parameter müssen die Bedingungen zur Unterscheidung des Regimes erfüllen. Wenn Sie beispielsweise zwischen risikoarmen und risikoreichen Regimen unterscheiden, müssen die vorgeschlagenen $ \ sigma_1 $ und $ \ sigma_2 $ $ \ sigma_1 <\ sigma_2 $ erfüllen. Beachten Sie, dass diese Bedingung für jedes einzelne Modell angegeben werden muss.

3. Implementierung des Markov-Switching-Modells durch Python

Implementieren Sie die Markov-Umschaltung für eine beliebige Anzahl von Regimen.

import numpy as np
import pandas as pd
from scipy import stats

import matplotlib.pyplot as plt
import japanize_matplotlib
%matplotlib inline

from tqdm.notebook import tqdm

regime = 3

#Anfangsparameter der Normalverteilung, der jedes Regime folgt
mu = [-1.4,-0.6,0]
sigma = [0.07,0.03,0.04]

Anstatt die Übergangswahrscheinlichkeit des Regimes selbst mit MCMC zu schätzen, wird die Übergangswahrscheinlichkeit als logistischer Funktionswert verwendet und ihre Argumente werden mit MCMC geschätzt. Zu diesem Zeitpunkt wird eine Logistikfunktion auf die nicht diagonale Komponente angewendet, und die diagonale Komponente erzeugt eine Übergangswahrscheinlichkeitsmatrix als Co-Ereignis. Behandeln Sie die Argumente der Logistikfunktion separat als "gen_prob" und den Wert der Logistikfunktion als "prob".

gen_prob = np.ones((regime,regime))*(-3)
gen_prob = gen_prob - np.diag(np.diag(gen_prob))
#Diagonale Komponente ist 0, nicht diagonale Komponente ist-Quadratische Matrix von 3

prob = np.exp(gen_prob)/(1+np.exp(gen_prob))
#Logistikfunktion anwenden

prob = prob + np.diag(1-np.dot(prob.T,np.ones((regime,1))).flatten())
#Als Restereignis der nicht diagonalen Komponente wird die Wahrscheinlichkeit der diagonalen Komponente erhalten und als Übergangswahrscheinlichkeitsmatrix verwendet.

#Eine Funktion zum Ermitteln der stationären Wahrscheinlichkeit der Markov-Kette
def cal_stationary_prob(prob, regime):
    # prob: 2-d array, shape = (regime, regime),Übergangswahrscheinlichkeitsmatrix
    # regime: int,Anzahl der Regime
    
    A = np.ones((regime+1, regime))
    A[:regime, :regime] = np.eye(regime)-prob
    
    return np.dot(np.linalg.inv(np.dot(A.T,A)),A.T)[:,-1]

#Funktion zur Berechnung der Protokollwahrscheinlichkeit
def cal_logL(x, mu, sigma, prob, regime):
    # x: 1-d array or pandas Series,Daten in chronologischer Reihenfolge
    # mu: 1=d array, len = regime,Anfangswert des Durchschnitts der Normalverteilung, der jedes Regime folgt
    # sigma: 1=d array, len = regime,Anfangswert der Varianz der Normalverteilung, der jedes Regime folgt
    # prob: 2-d array, shape = (regime, regime),Übergangswahrscheinlichkeitsmatrix
    # regime: int,Anzahl der Regime
    
    likelihood = stats.norm.pdf(x=x, loc=mu[0], scale=np.sqrt(sigma[0])).reshape(-1,1)
    for i in range(1,regime):
        likelihood = np.hstack([likelihood, stats.norm.pdf(x=x, loc=mu[i], scale=np.sqrt(sigma[i])).reshape(-1,1)])
    
    prior = cal_stationary_prob(prob, regime)
    
    logL = 0
    for i in range(length):
        temp = likelihood[i]*prior
        sum_temp = sum(temp)

        posterior = temp/sum_temp

        logL += np.log(sum_temp)
        
        prior = np.dot(prob, posterior)
        
    return logL

#Eine Funktion, die die Wahrscheinlichkeit berechnet, dass jeder Zeitpunkt zu jedem Regime gehört
def prob_regime(x, mu, sigma, prob, regime):
    # x: 1-d array or pandas Series,Daten in chronologischer Reihenfolge
    # mu: 1=d array, len = regime,Anfangswert des Durchschnitts der Normalverteilung, der jedes Regime folgt
    # sigma: 1=d array, len = regime,Anfangswert der Varianz der Normalverteilung, der jedes Regime folgt
    # prob: 2-d array, shape = (regime, regime),Übergangswahrscheinlichkeitsmatrix
    # regime: int,Anzahl der Regime
    
    likelihood = stats.norm.pdf(x=x, loc=mu[0], scale=np.sqrt(sigma[0])).reshape(-1,1)
    for i in range(1,regime):
        likelihood = np.hstack([likelihood, stats.norm.pdf(x=x, loc=mu[i], scale=np.sqrt(sigma[i])).reshape(-1,1)])
    
    prior = cal_stationary_prob(prob, regime)
    
    prob_list = []
    for i in range(length):
        temp = likelihood[i]*prior

        posterior = temp/sum(temp)
        
        prob_list.append(posterior)
        
        prior = np.dot(prob, posterior)
    
    return np.array(prob_list)

#Eine Funktion, die Parameter generiert, die von MCMC aktualisiert werden sollen
def create_next_theta(mu, sigma, gen_prob, epsilon, regime):
    # mu: 1=d array, len = regime,Anfangswert des Durchschnitts der Normalverteilung, der jedes Regime folgt
    # sigma: 1=d array, len = regime,Anfangswert der Varianz der Normalverteilung, der jedes Regime folgt
    # gen_prob: 2-d array, shape = (regime, regime),Argumente für logistische Funktionen
    # epsilon: float,Aktualisieren Sie die Breite des vorgeschlagenen Parameters
    # regime: int,Anzahl der Regime
    
    new_mu = mu.copy()
    new_sigma = sigma.copy()
    new_gen_prob = gen_prob.copy()
    
    new_mu += (2*np.random.rand(regime)-1)*epsilon
    new_mu = np.sort(new_mu)
    new_sigma = np.exp(np.log(new_sigma) + (2*np.random.rand(regime)-1)*epsilon)
    new_sigma = np.sort(new_sigma)[[2,0,1]]
    new_gen_prob += (2*np.random.rand(regime,regime)-1)*epsilon*0.1

    new_gen_prob = new_gen_prob - np.diag(np.diag(new_gen_prob))
    new_prob = np.exp(new_gen_prob)/(1+np.exp(new_gen_prob))
    new_prob = new_prob + np.diag(1-np.dot(new_prob.T,np.ones((regime,1))).flatten())
    
    return new_mu, new_sigma, new_gen_prob, new_prob

#Funktion, die MCMC ausführt
def mcmc(x, mu, sigma, gen_prob, prob, epsilon, trial, regime):
    # x: 1-d array or pandas Series,Daten in chronologischer Reihenfolge
    # mu: 1=d array, len = regime,Anfangswert des Durchschnitts der Normalverteilung, der jedes Regime folgt
    # sigma: 1=d array, len = regime,Anfangswert der Varianz der Normalverteilung, der jedes Regime folgt
    # gen_prob: 2-d array, shape = (regime, regime),Argumente für logistische Funktionen
    # epsilon: float,Aktualisieren Sie die Breite des vorgeschlagenen Parameters
    # trial: int,Anzahl der MCMC-Ausführungen
    # regime: int,Anzahl der Regime
    
    mu_list = []
    sigma_list = []
    prob_list = []
    logL_list = []
    mu_list.append(mu)
    sigma_list.append(sigma)
    prob_list.append(prob)
    
    for i in tqdm(range(trial)):
        new_mu, new_sigma, new_gen_prob, new_prob = create_next_theta(mu, sigma, gen_prob, epsilon, regime)
        
        logL = cal_logL(x, mu, sigma, prob, regime)
        next_logL = cal_logL(x, new_mu, new_sigma, new_prob, regime)
        
        ratio = np.exp(next_logL-logL)
        logL_list.append(logL)

        if ratio > 1:
            mu, sigma, gen_prob, prob = new_mu, new_sigma, new_gen_prob, new_prob

        elif ratio > np.random.rand():
            mu, sigma, gen_prob, prob = new_mu, new_sigma, new_gen_prob, new_prob

        mu_list.append(mu)
        sigma_list.append(sigma)
        prob_list.append(prob)
            
        if i%100==0:
            print(logL)
    
    return np.array(mu_list), np.array(sigma_list), np.array(prob_list), np.array(logL_list)

4. Datenbeispiel

TEPCOs überschüssige Gewinnspanne

Der Zeitraum reicht von 2003 bis 2019. Der Marktaktienkurs war topix, und die risikofreien Vermögenswerte waren japanische Staatsanleihen (10 Jahre).

Die Zeitreihendaten der Überschussgewinnrate sind Pandas 'DataFramex.

regime = 3

mu = [-1.4,-0.6,0]
sigma = [0.07,0.03,0.04]

gen_prob = np.ones((regime,regime))*(-3)
gen_prob = gen_prob - np.diag(np.diag(gen_prob))
prob = np.exp(gen_prob)/(1+np.exp(gen_prob))
prob = prob + np.diag(1-np.dot(prob.T,np.ones((regime,1))).flatten())

trial = 25000
epsilon = 0.1

mu_list, sigma_list, prob_list, logL_list = mcmc(x, mu, sigma, gen_prob, prob, epsilon, trial, regime)

prob_series = prob_regime(x, mu_list[-1], sigma_list[-1], prob_list[-1], regime)

Visualisieren Sie die Wahrscheinlichkeit jedes Regimes als gestapeltes Balkendiagramm.

fig = plt.figure(figsize=(10,5),dpi=200)
ax1 = fig.add_subplot(111)
ax1.plot(range(length), x, linewidth = 1.0, label="return")

ax2 = ax1.twinx()
for i in range(regime):
    ax2.bar(range(length), prob_series[:,i], width=1.0, alpha=0.4, label=f"regime{i+1}", bottom=prob_series[:,:i].sum(axis=1))

h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, bbox_to_anchor=(1.05, 1), loc='upper left')

ax1.set_xlabel('t')
ax1.set_ylabel(r'Übermäßige Rentabilität')
ax2.set_ylabel(r'Wahrscheinlichkeit')
ax1.set_title(f"Wahrscheinlichkeit des stark verteilten Regimes von TEPCO, epsilon={epsilon}, trial={trial}")

plt.subplots_adjust(left = 0.1, right = 0.8)
#plt.savefig("probability_regime2.png ")
plt.show()

In der folgenden Grafik ist das Regime für hohe Diversifikation und niedrige Rentabilität blau, das Regime für niedrige Diversifikation und niedrige Rentabilität gelb und das Regime für mittlere Diversifikation und hohe Rentabilität grün. probability_regime_with_excess_return.png

Die folgende Grafik zeigt den tatsächlichen Aktienkurs und nicht die Gewinnspanne. Es ist ersichtlich, dass der Zusammenbruch der Aktienkurse aufgrund des Atomunfalls nach dem Erdbeben im Großen Osten Japans eher ein Sprung als ein Regimewechsel ist. probability_regime_stock_price.png

Verweise

Recommended Posts

Markov-Schaltmodell von Python
Erklärung des Produktionsoptimierungsmodells durch Python
Python-Implementierung eines kontinuierlichen Hidden-Markov-Modells
Primzahlbeurteilung durch Python
Kommunikationsverarbeitung durch Python
Beamformer-Antwort von Python
Sprachvorhersagemodell von TensorFlow
Visualisieren Sie das Keras-Modell mit Python 3.5
Spracherkennung durch Python MFCC
EXE Web API von Python
Newcomer Trainingsprogramm von Python
Parametereinstellung durch Python Configparser
Pin Python von Conda verwaltet
Keyword-Extraktion mit MeCab (Python)
Pokemon-Klassifizierung nach Themenmodell
Zahlen durch 3 Ziffern trennen (Python)
Bildverarbeitung mit Python (Pillow)
Python wurde von C-Programmierern gestartet
Himbeere pi 1 Modell b, Python
Plattform (OS) Beurteilung durch Python
Sortieren nach Datum in Python
[Python] Sortierbar nach mehreren Bedingungen sortieren
Zusammenfassung des maschinellen Lernens von Python-Anfängern
Lerne Python durch Zeichnen (Turtle Graphics)
Python-Entwicklung unterstützt durch Jenkins-Unit-Test
Primzahlgenerator von Python
Python SQL-Anweisung Nach Zeit extrahieren
Betriebssystembestimmung durch Makefile mit Python
Geschätztes Probit-Modell nach binärem Antwortmodell
Modell generiert von Variational Autoencoder (VAE)
[Python] Gemischtes Gaußsches Modell mit Pyro
Abschnittsplanung Lernnotiz ~ von Python ~
Allgemeines Gaußsches Zustandsraummodell in Python
Verhalten von Python3 durch Sakuras Server
100 Sprachverarbeitung Knock Kapitel 1 von Python
Verwenden Sie das von fastText trainierte Modell von Python
Geschichte der Potenznäherung von Python
Benutzerdefiniertes Zustandsraummodell in Python
Sortieren von Dateien nach Namenskonvention mit Python