Hi all !
I have to create a regime switching model. I found this great model implemented by @junpenglao (here).
However, my dataset is larger and it takes too much time to run.
Here is the model:
with pm.Model() as m:
# Transition matrix
p = pm.Beta('p', alpha=10., beta=2., shape=2)
P = tt.diag(p)
P = tt.set_subtensor(P[0, 1], 1 - p[0])
P = tt.set_subtensor(P[1, 0], 1 - p[1])
# eta
alpha = pm.Normal('alpha', mu=0., sd=.1, shape=2)
sigma = pm.HalfCauchy('sigma', beta=1., shape=2)
eta1 = tt.exp(pm.Normal.dist(mu=alpha[0], sd=sigma[0]).logp(yshared))
y_tm1_init = pm.Normal('y_init', mu=0., sd=.1)
pNormal = pm.Bound(pm.Normal, lower=0.0)
rho = pNormal('rho', mu=1.0, sd=.1, testval=1.0)
eta2 = tt.zeros_like(eta1)
eta2 = tt.set_subtensor(eta2[0], tt.exp(
pm.Normal.dist(mu=alpha[1] + rho * y_tm1_init, sd=sigma[1]).logp(yshared[0])))
eta2 = tt.set_subtensor(eta2[1:], tt.exp(
pm.Normal.dist(mu=alpha[1] + rho * yshared[:-1], sd=sigma[1]).logp(yshared[1:])))
eta = tt.stack([eta1, eta2], axis=1)
xi_init = pm.Beta('xi_init', alpha=2., beta=2.)
ft_out = theano.shared(0.) # place holder
([ft, xi], updates) = theano.scan(ft_xit_dt,
sequences=eta,
outputs_info=[ft_out, xi_init],
non_sequences=P)
Xi = pm.Deterministic('Xi', xi)
# likelihood `target += sum(log(f))`
pm.Potential('likelihood', tt.sum(tt.log(ft)))
with m:
trace = pm.sample(5000, tune=5000, cores=3, init="adapt_diag")
I would like to break this long dataserie into a smaller list of fixed size ones, like a moving window.
So when I call the sampling function, it works on the whole list.
Furthermore, I would like to have for each series: it draws 20 times and keeps the one with the most likelyhood.
I haven’t used pymc3 much yet.
Should I create a custom walk function as shown here?
Thanks in advance for the time spent reading me and guiding me in the right direction.
Pierre