I’d like to do an adstock transformation (an econometric transformation used in advertising effectiveness modelling) within `pm.Deterministic()`.

There is lots of resources on how to do it (A Bayesian Approach to Media Mix Modeling by Michael Johns & Zhenyu Wang) but I couldn’t find anything on how to incorporate into a hierarchical model. Well, there is this article (https://towardsdatascience.com/bayesian-hierarchical-marketing-mix-modeling-in-pymc-684f6024e57a) but it does a for loop within `pm.Model()` which I tried but it makes it very slow (8 hours 32,000 draws).

My model looks like this:

``````coords = {"group": group_labels}
with pm.Model(coords=coords) as hierarchical:
# Hyperpriors
alpha_mu = pm.HalfNormal("alpha_mu", sigma=1)
beta_mu = pm.HalfNormal("beta_mu", sigma=1)
sigma_mu = pm.HalfNormal("sigma_mu", sigma=1)
...
# Priors
theta = pm.Beta("theta", a=alpha_mu, b=beta_mu, dims="groups")
b_spend = pm.HalfNormal("b_spend", sigma=sigma_mu, dims="groups")
...
# Data
group_idx = pm.MutableData("group", df["group_idx"], dims="obs_id")
spend = pm.MutableData("spend", df["spend"], dims="obs_id")
...
# Model
...
# Likelihood
...
``````

I took the `geometric_adstock` from here.

``````def geometric_adstock(x, alpha: float = 0.0, l_max: int = 12):
cycles = [
at.concatenate(
[at.zeros(i), x[: x.shape[0] - i]]
)
for i in range(l_max)
]
x_cycle = at.stack(cycles)
w = at.as_tensor_variable([at.power(alpha, i) for i in range(l_max)])
return at.dot(w, x_cycle)
``````

But of course it doesn’t work in my model.

I’m not an expert of vectorized operations nor linear algebra so this was my best attempt to make it work.

``````import aesara.tensor as at
import aesara

x_ = x.reshape((n_groups, -1))
X = at.stack(
[
at.stack(
[
at.concatenate([np.zeros(i), x_[r, : x_.shape[1] - i]])
for i in range(max_len)
]
)
for r in range(n_groups)
]
)
W = at.transpose(at.stack([at.power(theta, i) for i in range(max_len)]))
adstock_x = at.transpose(at.diagonal(at.dot(W, X), axis1=0, axis2=1))
``````

It’s ugly but it does work I think based on this toy data

``````# one theta per group
theta = np.array([
0.5,
0.25,
0.5,
0.1
])
# 4 obs per group
x = np.array([
0, 100, 0, 0, 0, 50, 0, 0, 0, 100, 0, 0, 100, 0, 0, 0
])
s = at.vector("s")
t = at.vector("t")
UPDATE: I’ve noticed that because of my implementation of `at_adstock`, `pm.sample()` automatically switches to Metropolis for “theta” and the resulting trace is non-sensical. When I force the step to be NUTS I get an error:
`NullTypeGradError: `grad` encountered a NaN. This variable is Null because the grad method for input 0 (Reshape{3}.0) of the ExtractDiag{offset=0, axis1=0, axis2=1, view=False} op is not implemented.` BTW I’ve replaced `at.diagonal` with `pymc.math.extract_diag` hoping it would fix it. It didn’t. I know this is all due to my dirty implementation of `at_adstock`. I just can’t think of a better way of doing it.