我正在尝试按照本示例的思路在 pymc3 中实现 GARCH 模型。为此,我尝试如下实现 GARCH(1, 1) 分布
import pymc3 as pm
from pymc3.distributions import Continuous, Normal
class GARCH(Continuous):
def __init__(self, alpha_0=None, alpha_1=None, beta_1=None, sigma_0=None, *args, **kwargs):
super(GARCH, self).__init__(*args, **kwargs)
self.alpha_0 = alpha_0
self.alpha_1 = alpha_1
self.beta_1 = beta_1
self.sigma_0 = sigma_0
self.mean = 0
def logp(self, values):
sigma = self.sigma_0
alpha_0 = self.alpha_0
alpha_1 = self.alpha_1
beta_1 = self.beta_1
x_prev = values[0]
_logp = Normal.dist(0., sd=sigma).logp(x_prev)
for x in values[1:]:
sigma = pm.sqrt(alpha_0 + alpha_1 * (x_prev/sigma)**2
+ beta_1 * sigma**2)
_logp = _logp + pm.Normal(0., sd=sigma).logp(x)
x_prev = x
return _logp
澄清一下,这是 GARCH(1,1) 模型的对数似然。波动过程是一个时间序列,其中时间 t 的波动取决于时间 t-1 的残差。但是要确定时间 t-1 的残差,我们需要时间 t-1 的波动率。
无论如何,这对我的问题并不重要。重要的是不能通过矢量化 for 循环来计算可能性(这是在帖子顶部的链接中完成的)。因此,您需要一个显式循环,它在每一步首先更新波动率,然后确定观察到的回报的可能性。
但是上面的代码不起作用。如果我尝试建立一个像
import numpy as np
returns = np.genfromtxt("SP500.csv")[-200:]
garchmodel = pm.Model()
with garchmodel:
alpha_0 = pm.Exponential('alpha_0', 30., testval=.02)
alpha_1 = pm.Uniform('alpha_1', lower=0, upper=1, testval=.9)
upper = pm.Deterministic('upper', 1-alpha_1)
beta_1 = pm.Uniform('beta_1', lower=0, upper=upper, testval=.05)
sigma_0 = pm.Exponential('sigma_0', 30., testval=.02)
garch = GARCH('garch', alpha_0=alpha_0, alpha_1=alpha_1,
beta_1=beta_1, sigma_0=sigma_0, observed=returns)
“SP500.csv”文件可以在例如github上找到
此代码生成错误:
ValueError: length not known
我很确定这是因为 for 循环与 theano 冲突。我该如何处理?