Calibración del modelo de precio del petróleo con Filtro de Kalman y MLE en python

Question

Estoy intentando calibrar un proceso de reverberación media de un factor en Python 3. El proceso está definido como:

\begin{equation} dX = k(\alpha - X)dt + \sigma dW , \end{equation}

donde $\alpha = \mu - \frac{\sigma^2}{2k}$ es el precio logarítmico a largo plazo y $k$ es la velocidad de ajuste. Bajo la probabilidad neutral al riesgo $Q$ escribimos la ecuación anterior como

\begin{equation} dX = k(\alpha^* - X)dt + \sigma dW^*. \end{equation}

Aquí $\alpha^*= \alpha - \lambda$, donde $\lambda$ es el precio del riesgo de mercado. Basado en el Teorema de Girsanov, $dW^*$ es una Movimiento Browniano bajo la medida martingala $Q$.

Mis ecuaciones de medición y transición son:

• La ecuación de medición relaciona la serie de tiempo de variables observables, en mi caso los precios futuros para diferentes vencimientos, con la variable de estado no observable, el precio spot:

\begin{equation} y_{t} = Z_{t}X_{t} + d_{t} + \epsilon_{t}, \qquad t = 1, ..., NT \end{equation}

donde

$y_{t}=\left[\ln F\left(T_{i}\right)\right]$, $i=1,...,N$, $N\times 1$ vector de observables,

$d_{t}=\left[\left(1-e^{-\kappa T_{1}}\right) \alpha^{*}+\frac{\sigma^{2}}{4 \kappa}\left(1-e^{-2 k T_{1}}\right)\right], \quad i=1, \ldots, N, \quad N \times 1$ vector,

$Z_{t}=\left[e^{-\kappa T_{i}}\right], \quad i=1, \ldots, N, \quad N \times 1$ vector,

$\epsilon_{t}, \quad N\times 1$ vector de perturbaciones no correlacionadas en serie con $\mathbb{E}(\epsilon_{t}) = 0$ y $Var(\epsilon_{t}) = H$.

• La ecuación de transición es una versión discreta en tiempo del proceso estocástico de precios del petróleo O-U:

\begin{equation} X_{t}=c_{t}+Q_{t} X_{t-1}+\eta_{t}, \quad t=1, \ldots, N T \label{eq:38}, \end{equation}

donde

$c_{t}=\kappa \alpha \Delta t$

$Q_{t}=1-\kappa \Delta t$

$\eta_{t}$, perturbaciones no correlacionadas en serie con $\mathbb{E}(\eta_{t}) = 0$ y $Var(\eta_{t}) = \sigma^2 \Delta t$.

El objetivo es calibrar los parámetros del modelo utilizando la serie de tiempo de precios futuros de diferentes vencimientos a través de MLE. ¡Cualquier ayuda con la parte de programación (en Python 3) sería muy apreciada!

Answer 1

Ampliando la respuesta de @ir7, aquí hay algún código/pseudocódigo de pykalman para ayudarte a empezar. Esto se puede ajustar de muchas maneras, pero he dejado algunos parámetros para darte una idea. También dejé un enlace a la documentación al final. Las funciones configurarán Filtros de Kalman que se aplican a tus datos y posteriormente esos datos se alimentan a una regresión que posteriormente se llama en una función de backtest:

import pandas as pd
import numpy as np
from pykalman import KalmanFilter

def your_function_name(x):
    your_filter_name = KalmanFilter(transition_matrices = [1],
                                    observation_matrices = [1],
                                    observation_covariance=1,
                                    transition_covariance=.01,
                                    initial_state_mean = 0,
                                    initial_state_covariance = 1)

    state_means, _ = kf.filter(x.values)
    state_means = pd.Series(state_means.flatten(), index=x.index)
    return state_means

def your_regression_filter(x, y):
    delta = 1e-3
    trans_cov = delta / (1 - delta) * np.eye(2)  #random walk wiggle
    obs_mat = np.expand_dims(np.vstack([[x], [np.ones(len(x))]]).T, axis=1)
    kf = KalmanFilter(n_dim_obs=1,
                      n_dim_state=2,
                      initial_state_mean=[0,0],
                      initial_state_covariance=np.ones((2, 2)),
                      transition_matrices=np.eye(2),
                      observation_matrices=obs_mat,
                      observation_covariance=2,
                      transition_covariance=trans_cov)

    # Usa las observaciones y para obtener estimaciones en ejecución y errores para los parámetros de estado
    state_means, state_covs = kf.filter(y.values)
    return state_means

Con estas 2 funciones definirías una función de backtest en la que extraerías state_means llamando:

state_means = your_regression_filter(your_function_name(x), your_function_name(y))

Cómo uses state_means a partir de aquí depende de ti.

Para más información: documentación de pykalman

Answer 2

Un recurso que tiene Filtro y Suavizador de Kalman, así como algoritmos de Expectation-Maximization para un Modelo Gaussiano Lineal es el módulo pykalman. También puedes revisar el módulo statsmodels.