Оценка модели AR с переменным коэффициентом с помощью statsmodels - PullRequest
Новая
       86

Оценка модели AR с переменным коэффициентом с помощью statsmodels

0 голосов
/ 29 апреля 2020

Я вхожу в модели пространства состояний и их оценку в пакете python statsmodels. Я получил пример из руководства пользователя для работы, но не смог получить разумные результаты для своего собственного примера, и я не совсем понимаю, почему.

Я хочу оценить AR (1) с изменяющимися во времени коэффициентами, то есть y_t = a_t + b_t*y_{t-1} + e_t (уравнение измерения), где вектор коэффициентов alpha_t := (a_t,b_t)' следует за VAR (1) (уравнение перехода).

Здесь симуляция 

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

# simulate states
T = 500 # nobs
A = np.diag([0.8,0.8]) # transition matrix
c = np.array([4,0.1]) # state intercept
Q = np.eye(2)
Q[1,1] /= 1000 # cov of state error

# draw state variable errors
v = np.random.multivariate_normal([0,0],Q,size=(T,))
init = np.array([20,0.5])
alpha = np.vstack([init[None,:],np.zeros((T,2))])
for i in range(T):
    alpha[i+1,:] = c + A@alpha[i,:] + v[i,:]

# plot processes
fig, ax = plt.subplots(2,2,figsize=(17,8))
for i in range(2):
    ax[i,0].plot(alpha[1:,i], label=r'$\alpha_{i}$'.format(i=i+1))
    ax[i,0].legend(loc=1)
    ax[i,1].plot(v[1:,i], label=r'$v_{i}$'.format(i=i+1))
    ax[i,1].legend(loc=1)

plt.show()

# simulate y_t
e = np.random.randn(T)*5
y = np.append(40,np.zeros(T))
for i in range(T):
    y[i+1] = alpha[i+1,0] + alpha[i+1,1]*y[i] + e[i]

# plot y_t
fig, ax = plt.subplots(1,1,figsize=(12,5))
ax.plot(y)
ax.set_title(r'$y_t$')
plt.show()

Я пытался выбрать параметры так, чтобы y_t вел себя стабильно, но я думаю, что это не гарантировано, поэтому лучше составьте график процессов перед оценкой.

Вот как я представляю модель в statsmodels: 

class ARVarCoefModel(sm.tsa.statespace.MLEModel):

    param_names = ['sigma2.e','c0','c1','a11','a22','sigma2.v1','sigma2.v2'] # for output readability
    start_params = [25,4,0.1,0.8,0.8,1,0.001] # initiate estimation at true parameters

    def __init__(self, endog, design_mat):
        k_states = k_posdef = design_mat.shape[1]
        super().__init__(endog, k_states=k_states, k_posdef=k_posdef,
            initialization='approximate_diffuse',loglikelihood_burn=k_states) #approximate_diffuse

        self['design'] = exog.T[None,:,:] 
        self['selection'] = np.eye(2) 

        self.diag_idx = np.diag_indices(k_states)


    # variance must be positive; keep transition matrix params between (0,1) 
    def transform_params(self,params):
        constrained = params.copy()
        constrained[0] = constrained[0]**2 
        constrained[3] = 1/(1+np.exp(constrained[3])) 
        constrained[4] = 1/(1+np.exp(constrained[4])) 
        constrained[5] = constrained[5]**2 
        constrained[6] = constrained[6]**2 
        return constrained 

    def untransform_params(self,params):
        unconstrained = params.copy()
        unconstrained[0] = unconstrained[0]**0.5 
        unconstrained[3] = np.log(1/unconstrained[3] -1) 
        unconstrained[4] = np.log(1/unconstrained[4] -1) 
        unconstrained[5] = unconstrained[5]**0.5 
        unconstrained[6] = unconstrained[6]**0.5
        return unconstrained 

    def update(self,params, **kwargs):
        params = super().update(params, **kwargs)

        self['obs_cov',0,0] = params[0]
        self['state_intercept'] = params[1:3]
        self[('transition',) + self.diag_idx] = params[3:5]
        self[('state_cov',) + self.diag_idx] = params[5:]

Теперь оцените 

endog = y[1:]
exog = np.hstack([np.ones((len(endog),1)),y[:-1,None]])
model = ARVarCoefModel(endog,exog)
model.fit(method='bfgs',maxiter=1000)

Независимо от того, как часто я запускаю симуляцию, какой метод оптимизации я использую, оптимизация сходится, но результаты имеют мало смысла. По какой-то причине (или, может быть, это совпадение) параметры a11 и sigma2.v11 резко занижены (относительно близки к нулю).

Сглаженные состояния также не приближаются к тому, какими они должны быть: 

fig, ax = plt.subplots(2,1,figsize=(14,7))
for i in range(2):

    ax[i].plot(alpha[1:,i],label=r'$\alpha_{}$'.format(i+1))
    ax[i].plot(res.smoothed_state[i,1:],'g--',label=r'$\hat \alpha_{}$ (smoothed)'.format(i+1))

    ax[i].legend()

plt.show()

У меня нет практического опыта работы с моделями пространства состояний в этой более общей форме, поэтому, возможно, я выбрал плохой пример и это просто некорректная задача. Но я бы подумал, что это должен быть довольно стандартный пример (подобный практический пример я видел в Lutkepohl (2005) «Новое введение в анализ множественных временных рядов», глава 18.5, которую мне не удалось воспроизвести). Поэтому я надеюсь, что просто делаю что-то не так, представляя модель в statsmodels. (Кстати, это версия 0.9.0; Python 3.6.5)

...