Modèle de commutation de Markov par Python

1. "Régime" des données économiques

Lorsqu'un changement soudain et important survient dans l'économie, tel que l'effondrement de la bulle informatique ou le choc Lehman, on peut considérer qu'un état a fondamentalement changé en arrière-plan. Un tel modèle de changement d'état fondamental est appelé «modèle de changement de régime». La modélisation d'un changement de régime peut être difficile si l'état modifié n'est pas observable. Cependant, la modélisation est possible en donnant la distribution de probabilité du régime avec des paramètres et en traitant la transition entre les régimes comme une chaîne de Markov. Les paramètres qui maximisent la vraisemblance (moyenne / écart type de la distribution normale, probabilité de transition du régime) peuvent être obtenus par MCMC.

2. Explication du modèle de commutation de Markov

Régression linéaire avec régime

Considérons le modèle de changement de régime du modèle de régression linéaire suivant.

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

Les régimes $ 1 $ à $ n $ à chaque instant $ t $ sont exprimés comme suit.

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

Si le régime peut être observé, l'équation de régression s'exprime comme suit avec le régime en indice.

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

A ce moment, $ y \ left (t \ right) $ est en moyenne $ x \ left (t \ right) \ beta_ {I \ left (t \ right)} $, distribué $ \ sigma ^ 2_ {I \ left ( Il suit une distribution normale de t \ right)} $. Appliquez la méthode la plus probable qui maximise la vraisemblance logarithmique. Autrement dit, le paramètre qui maximise la fonction de vraisemblance suivante est obtenu.

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

Si le régime ne peut pas être observé directement, alors $ I \ left (t \ right) $ est également considéré comme suivant une distribution de probabilité conditionnelle dans l'histoire passée. Soit $ \ mathcal {F} \ left (t \ right) $ l'information obtenue au temps $ t $. Notez que $ f \ left (t \ right) $ n'est pas conditionné sur $ \ mathcal {F} \ left (t \ right) $.

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

Par conséquent, la probabilité logarithmique à maximiser est la suivante.

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

Chaîne de Markov

Distribution de probabilité de $ y \ left (t \ right) $ en donnant $ f \ left (I \ left (t \ right) | \ mathcal {F} \ left (t-1 \ right) \ right) $ Est décidé. Supposons que $ I \ left (t \ right) $ satisfait la propriété Markov du premier ordre.

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

Probabilité de transition de l'état j à l'état i $ P \ left (I \ left (t \ right) = i | I \ left (t-1 \ right) = j \ right) $ à $ p_ {i, j} $ Et considérez la matrice de probabilité de transition suivante.

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

Probabilité conditionnelle du régime $ P \ left (I \ left (t \ right) = j | \ mathcal {F} \ left (t \ right) \ right) étant donné les informations jusqu'au temps $ t $ $ Est calculé par mise à jour bayésienne.

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

La probabilité initiale est la probabilité stationnaire de la chaîne de Markov. On sait que la matrice de probabilité de transition $ P $, la matrice unitaire $ I $ et le vecteur $ 1 $ avec tous les éléments de $ 1 $ sont utilisés, et la probabilité stationnaire $ \ pi ^ {\ ast} $ peut être obtenue comme suit.

\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ère rangée
\end{eqnarray}

Chaîne de Markov Méthode de Monte Carlo

La méthode Markov Chain Monte Carlo (MCMC) est utilisée comme méthode pour obtenir les paramètres qui maximisent la vraisemblance logarithmique. L'explication est omise ici.

Les paramètres proposés doivent répondre aux conditions de distinction du régime. Par exemple, lors de la distinction entre les régimes à faible risque et à haut risque, les $ \ sigma_1 $ et $ \ sigma_2 $ proposés doivent satisfaire $ \ sigma_1 <\ sigma_2 $. Notez que cette condition doit être spécifiée pour chaque modèle individuel.

3. Implémentation du modèle de commutation Markov par Python

Implémentez la commutation de Markov pour n'importe quel nombre de régimes.

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

#Paramètres initiaux de la distribution normale que chaque régime suit
mu = [-1.4,-0.6,0]
sigma = [0.07,0.03,0.04]

Au lieu d'estimer la probabilité de transition du régime lui-même avec MCMC, la probabilité de transition est utilisée comme valeur de fonction logistique et ses arguments sont estimés avec MCMC. A ce moment, une fonction logistique est appliquée à la composante hors diagonale, et la composante diagonale génère une matrice de probabilité de transition en tant que co-événement. Traitez les arguments de la fonction logistique comme «gen_prob» et la valeur de la fonction logistique comme «prob» séparément.

gen_prob = np.ones((regime,regime))*(-3)
gen_prob = gen_prob - np.diag(np.diag(gen_prob))
#La composante diagonale est 0, la composante hors diagonale est-Matrice carrée de 3

prob = np.exp(gen_prob)/(1+np.exp(gen_prob))
#Appliquer la fonction logistique

prob = prob + np.diag(1-np.dot(prob.T,np.ones((regime,1))).flatten())
#En tant qu'événement résiduel de la composante hors diagonale, la probabilité de la composante diagonale est obtenue et utilisée comme matrice de probabilité de transition.

#Une fonction pour trouver la probabilité stationnaire de la chaîne de Markov
def cal_stationary_prob(prob, regime):
    # prob: 2-d array, shape = (regime, regime),Matrice de probabilité de transition
    # regime: int,Nombre de régimes
    
    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]

#Fonction pour calculer la vraisemblance logarithmique
def cal_logL(x, mu, sigma, prob, regime):
    # x: 1-d array or pandas Series,Données par ordre chronologique
    # mu: 1=d array, len = regime,Valeur initiale de la moyenne de la distribution normale que suit chaque régime
    # sigma: 1=d array, len = regime,Valeur initiale de la variance de la distribution normale que chaque régime suit
    # prob: 2-d array, shape = (regime, regime),Matrice de probabilité de transition
    # regime: int,Nombre de régimes
    
    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

#Une fonction qui calcule la probabilité que chaque point temporel appartienne à chaque régime
def prob_regime(x, mu, sigma, prob, regime):
    # x: 1-d array or pandas Series,Données par ordre chronologique
    # mu: 1=d array, len = regime,Valeur initiale de la moyenne de la distribution normale que suit chaque régime
    # sigma: 1=d array, len = regime,Valeur initiale de la variance de la distribution normale que chaque régime suit
    # prob: 2-d array, shape = (regime, regime),Matrice de probabilité de transition
    # regime: int,Nombre de régimes
    
    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)

#Une fonction qui génère des paramètres à mettre à jour par MCMC
def create_next_theta(mu, sigma, gen_prob, epsilon, regime):
    # mu: 1=d array, len = regime,Valeur initiale de la moyenne de la distribution normale que chaque régime suit
    # sigma: 1=d array, len = regime,Valeur initiale de la variance de la distribution normale que chaque régime suit
    # gen_prob: 2-d array, shape = (regime, regime),Arguments de la fonction logistique
    # epsilon: float,Mettre à jour la largeur du paramètre proposé
    # regime: int,Nombre de régimes
    
    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

#Fonction qui exécute MCMC
def mcmc(x, mu, sigma, gen_prob, prob, epsilon, trial, regime):
    # x: 1-d array or pandas Series,Données par ordre chronologique
    # mu: 1=d array, len = regime,Valeur initiale de la moyenne de la distribution normale que chaque régime suit
    # sigma: 1=d array, len = regime,Valeur initiale de la variance de la distribution normale que chaque régime suit
    # gen_prob: 2-d array, shape = (regime, regime),Arguments de la fonction logistique
    # epsilon: float,Mettre à jour la largeur du paramètre proposé
    # trial: int,Nombre d'exécutions MCMC
    # regime: int,Nombre de régimes
    
    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. Exemple de données

Marge bénéficiaire excédentaire de TEPCO

La période est de 2003 à 2019. Le cours boursier était topix et les actifs sans risque étaient des obligations d'État japonaises (10 ans).

Les données de la série chronologique du taux de profit excédentaire sont 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)

Visualisez la probabilité de chaque régime sous forme de graphique à barres empilées.

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'Rentabilité excessive')
ax2.set_ylabel(r'probabilité')
ax1.set_title(f"Probabilité du régime hautement distribué de TEPCO, epsilon={epsilon}, trial={trial}")

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

Dans le graphique ci-dessous, le régime de diversification élevée et de faible rentabilité est bleu, le régime de faible diversification et de faible rentabilité est jaune et le régime de diversification moyenne et de rentabilité élevée est vert.

Le graphique ci-dessous trace le prix réel de l'action, pas la marge bénéficiaire. On peut voir que l'effondrement des cours des actions dû à l'accident nucléaire après le grand tremblement de terre au Japon oriental est présumé être un saut plutôt qu'un changement de régime.

Les références

--Takahiro Komatsu (2018), * Stratégie d'investissement optimale: théorie et pratique de la technologie de portefeuille *, Tokyo: Asakura Shoten --Tatsuyoshi Okimoto (2014), Application du modèle de changement de Markov à la macroéconomie et à la finance, * Journal of Japan Statistical Society *, 44 (1), 137-157