Python ARMA MLE (Реализация алгоритмов из литературы)
обзор
Я пытаюсь реализовать оптимизацию параметров авторегрессионного скользящего среднего (ARMA) с использованием оценки максимального правдоподобия (MLE) с помощью фильтра Калмана. Я знаю, что могу соответствовать моделям ARMA, используя пакет statsmodels в Python, но я хочу написать свою собственную реализацию вероятности ARMA и последующей оптимизации в качестве прототипа для будущей реализации C/C++. Кроме того, когда я просматриваю документацию по statsmodels, я обнаруживаю, что statsmodels Kalman Filter Log правдоподобия реализует несколько иное выражение, чем я нашел в литературе.
Алгоритмы
Чтобы вычислить вероятность регистрации ARMA, я следую статье Перлмана 1980 года:
Pearlman, JG "Алгоритм для точной вероятности процесса авторегрессии-скользящего среднего высокого порядка". Биометрика 67,1 (1980): 232-233.). Доступно от JSTOR.
Чтобы вычислить начальную матрицу P, я следую алгоритму в
Gardner G., Andrew C. Harvey и Garry DA Phillips. "Алгоритм AS 154: алгоритм точной оценки максимального правдоподобия моделей авторегрессии-скользящего среднего с помощью фильтрации Калмана". Журнал Королевского статистического общества. Серия C (Прикладная статистика) 29,3 (1980): 311-322. Доступно от JSTOR.
Для начальных значений параметров в настоящее время я использую внутренний метод, который модели statsmodels ARMA используют для вычисления первоначального предположения для параметров ARMA. В будущем я планирую перейти к собственной реализации, но я использую _fit_starts_params, пока отлаживаю MLE.
Для оптимизации MLE я просто использую решатель L-BFGS в Scipy.
Код
import numpy as np
import statsmodels.api as sm
import statsmodels.tsa.arima_model
import scipy.optimize
class ARMA(object):
def __init__(self, endo, nar, nma):
np.random.seed(0)
# endogenous variables
self.endo = endo
# Number of AR terms
self.nar = nar
# Number of MA terms
self.nma = nma
# "Dimension" of the ARMA fit
self.dim = max(nar, nma+1)
# Current ARMA parameters
self.params = np.zeros(self.nar+self.nma, dtype='float')
def __g(self, ma_params):
'''
Build MA parameter vector
'''
g = np.zeros(self.dim, dtype='float')
g[0] = 1.0
if self.nma > 0:
g[1:self.nma+1] = ma_params
return g
def __F(self, ar_params):
'''
Build AR parameter matrix
'''
F = np.zeros((self.dim, self.dim), dtype='float')
F[:self.nar, 0] = ar_params
for i in xrange(1, self.dim):
F[i-1, i] = 1.0
return F
def __initial_P(self, R, T):
'''
Solve for initial P matrix
Solves P = TPT' + RR'
'''
v = np.zeros(self.dim*self.dim, dtype='float')
for i in xrange(self.dim):
for j in xrange(self.dim):
v[i+j*self.dim] = R[i]*R[j]
R = np.array([R])
S = np.identity(self.dim**2, dtype='float')-np.kron(T, T)
V = np.outer(R, R).ravel('F')
Pmat = np.linalg.solve(S,V).reshape(self.dim, self.dim, order='F')
return Pmat
def __likelihood(self, params):
'''
Compute log likehood for a parameter vector
Implements the Pearlman 1980 algorithm
'''
# these checks are pilfered from statsmodels
if self.nar > 0 and not np.all(np.abs(np.roots(np.r_[1, -params[:self.nar]]))<1):
print 'AR coefficients are not stationary'
if self.nma > 0 and not np.all(np.abs(np.roots(np.r_[1, -params[-self.nma:]]))<1):
print 'MA coefficients are not stationary'
ar_params = params[:self.nar]
ma_params = params[-self.nma:]
g = self.__g(ma_params)
F = self.__F(ar_params)
w = self.endo
P = self.__initial_P(g, F)
n = len(w)
z = np.zeros(self.dim, dtype='float')
R = np.zeros(n, dtype='float')
a = np.zeros(n, dtype='float')
K = np.dot(F, P[:, 0])
L = K.copy()
R[0] = P[0, 0]
for i in xrange(1, n):
a[i-1] = w[i-1] - z[0]
z = np.dot(F, z) + K*(a[i-1]/R[i-1])
Kupdate = -(L[0]/R[i-1])*np.dot(F, L)
Rupdate = -L[0]*L[0]/R[i-1]
P -= np.outer(L, L)/R[i-1]
L = np.dot(F, L) - (L[0]/R[i-1])*K
K += Kupdate
R[i] = R[i-1] + Rupdate
if np.abs(R[i] - 1.0) < 1e-9:
R[i:] = 1.0
break
for j in xrange(i, n):
a[j] = w[j] - z[0]
z = np.dot(F, z) + K*(a[i-1]/R[i-1])
likelihood = 0.0
for i in xrange(n):
likelihood += np.log(R[i])
likelihood *= -0.5
ssum = 0.0
for i in xrange(n):
ssum += a[i]*a[i]/R[i]
likelihood += -0.5*n*np.log(ssum)
return likelihood
def fit(self):
'''
Fit the ARMA model by minimizing the loglikehood
Uses scipy.optimize.minimize
'''
sm_arma = statsmodels.tsa.arima_model.ARMA(endog=self.endo, order=(self.nar, self.nma, 0))
params = statsmodels.tsa.arima_model.ARMA._fit_start_params_hr(sm_arma, order=(self.nar, self.nma, 0))
opt = scipy.optimize.minimize(fun=self.__likelihood, x0=params, method='L-BFGS-B')
print opt
# Test the code on statsmodels sunspots data
nar = 2
nma = 1
endo = sm.datasets.sunspots.load_pandas().data['SUNACTIVITY'].tolist()
arma = ARMA(endo=endo, nar=nar, nma=nma)
arma.fit()
вопросы
Я считаю, что приведенный выше пример не сходится. При третьем вызове ARMA._likelihood код выдает следующее предупреждение:
RuntimeWarning: invalid value encountered in log
likelihood += np.log(R[i])
что происходит, потому что ARMA._initial_P решает для матрицы, где P[0][0] < 0.0. В этот момент текущие оценки параметров AR становятся нестационарными. Все последующие итерации затем предупреждают, что параметры AR и MA нестационарны.
Вопросы
Это правильная реализация? Я проверил, что исходная матрица P удовлетворяет уравнению, которому она должна удовлетворять. Для расчета вероятности я вижу несколько вариантов поведения, которые я ожидаю от статьи Перлмана:
- R стремится к одному. Для процесса чистого AR с параметрами p AR он достигает этого предела за p шагов. По сути, оператор break в функции _likelihood вступает в силу после итераций шагов алгоритма Перлмана.
- L стремится к нулевому вектору.
- K имеет тенденцию к Fg. Я проверяю это, глядя на abs(K - Fg) при расчете вероятности.
- После предупреждения об отрицательном значении в логарифме вышеуказанные пределы больше не соблюдаются.
Я также попытался реализовать преобразование параметров ARMA, чтобы предотвратить переполнение / недополнение, как рекомендовано в
Джонс, Ричард Х. "Подбор максимального правдоподобия моделей ARMA для временных рядов с отсутствующими наблюдениями". Technometrics 22.3 (1980): 389-395. Доступно от JSTOR.
Это преобразование, казалось, не влияло на ошибки, которые я наблюдал.
Если реализация верна, то как мне обработать отрицательные значения R? Кажется, проблема возникает, когда scipy.optimize возвращает вектор параметров, который соответствует матрице P, для которой верхний диагональный элемент является отрицательным. Предполагается, что подпрограмма оптимизации ограничена для предотвращения отрицательных значений R? Я также попытался использовать сложные логарифмы для отрицательных значений, а также изменить все числовые параметры dtype на 'complex'. Например:
def complex_log(val): ''' Complex logarithm for negative values Returns log(val) + I*pi ''' if val < 0.0: return complex(np.log(np.abs(val)), np.pi) return np.log(val)
Однако scipy.optimize не может обрабатывать комплексные функции, поэтому это предполагаемое исправление пока не работает. Любые рекомендации по предотвращению или обращению с этим поведением?
Спасибо за чтение этого далеко. Буду признателен за любую оказанную помощь.