The current model I am using is:
with pm.Model() as model:
alpha = pm.Normal('alpha', mu = 1, sd = 10)
BoundedDelta = pm.Bound(pm.Normal, lower = -1, upper = 1)
delta = BoundedDelta('delta', mu = 0, sd = 3)
var_beta = pm.InverseGamma('var_beta',alpha=.001, beta=.001)
var_return = pm.InverseGamma('var_return', alpha=.001, beta=.001)
intercept = pm.Normal('intercept', mu = 0, sd = 1)
betas = StochasticBeta('betas', mu=0, alpha=alpha, delta=delta, var=var_beta, shape=len(returns_data_copy['excess_sym']))
r_mu = T.prod(betas, returns_data_copy['excess_market'].values) + T.sum(intercept)
observed = pm.Normal('observed', mu=r_mu, tau=1/var_return, observed=returns_data_copy['excess_sym'].values)
trace = pm.sample(2000, tune=1000, cores=4, progressbar=True)
pm.summary(trace)
with StochasticBeta( ) as:
class StochasticBeta(pm.distributions.distribution.Continuous):
def __init__(self, var, mu, alpha, delta, init=pm.Flat.dist(), *args, **kwargs):
super(StochasticBeta, self).__init__(*args, **kwargs)
self.var = var = T.as_tensor_variable(var)
self.mu = mu = T.as_tensor_variable(mu)
self.alpha = alpha = T.as_tensor_variable(alpha)
self.delta = delta = T.as_tensor_variable(delta)
self.init = init
self.mode = T.as_tensor_variable(0.)
def logp(self, x):
var = self.var
mu = self.mu
alpha = self.alpha
delta = self.delta
init = self.init
var = self.var
x_im1 = x[:-1]
x_i = x[1:]
innov_like = pm.Normal.dist(mu=alpha + delta * (x_im1 - alpha) + mu, tau=1/var).logp(x_i)
return init.logp(x[0]) + T.sum(innov_like)
Are there any recommendations on optimizing performance? How would I consider a reparameterization or scaling? It’s averaging roughly 10mins to complete but there are frequent divergences in the chains or failure to converge at all.