Time-varying coefficients

[oliver-kutis]

Is there a plan for varying coefficients implementation? If so, do you have a specific timeline?
I don't expect that the implementation is possible to be similar to Orbit so I've tried to implement those on my own using pm.GaussianRandomWalk (or pm.MvGaussianRandomWalk). However, as I have been trying to do so, the model sampling becomes very inefficient or doesn't converge even if I let it run for a very long time with high number of samples - when trying to generate coefficients for 2 or more variables. I am still trying to implement this but it seems impossible with my current knowledge of PyMC.

If you have plans for this feature, would it be possible to share your thoughts on how to do this efficiently in PyMC? Thank you!

[ricardoV94]

Thanks for reaching out @Oliver131313

Others will be more apt to reply than me. I just wanted to ask if you can share more details of what you have tried?
That may help providing specific advice. Maybe there's an obvious improvement that could be done to your approach.

[twiecki]

Yes, we plan to add time-varying parameters using GPs. Time-line isn't set yet for that.

[juanitorduz]

I suspect you have already seen this blog post where we do time-varying coefficients with Gaussian random walks. I have not really experimented with more than 2 variables, but I suspect GPs are the better choice given some success stories from PyMC Labs. We definitely want to do it!

A good starting point would be to modify the example https://www.pymc-marketing.io/en/stable/notebooks/mmm/mmm_example.html so that in the data generating process we add time-varying coefficients and then you could create a PR with your current approach and we could see if there is a way to make it work :)

[oliver-kutis]

Thank you all for feedback. @juanitorduz - yes I have seen your post. To be more precise, I have denifed the model succesfully using MvGaussianRandomWalk class, however it has taken too long to sample from, that I've rather deleted it and didn't plan to use it.

Currently, I'm trying to model the contributions as GPs, however I can't find resources with implementations which would include doing this with multiple dimensions (as in my applications for real-word problems). I am quite new to bayesian modeling and pymc as well. I'm planning to learn a lot but it's hard to find any consistent materials (I prefer a course).

I've tried to modify the code from the library using GPs but I'm either:

• incorrectly specifying the model, or
• finding myself waiting a very long time for the sampling to finish

I'm really struggling :)

Here is some code:

`
date_var = 'date_week'
target_var = 'y'
delay_vars = ['x1', 'x2']
control_vars = ['event_1', 'event_2']
fourier_vars = ['sin_order_1', 'cos_order_1', 'sin_order_2', 'cos_order_2']

with pm.Model() as gp_model:

# ------------------ Coords ------------------
gp_model.add_coord('date', data_df[date_var].values, mutable=True)
gp_model.add_coord('delay_vars', delay_vars, mutable=True)
gp_model.add_coord('control_vars', control_vars, mutable=True)
gp_model.add_coord('fourier_vars', fourier_vars, mutable=True)

# ------------------ Data ------------------
delay_vars_data_ = pm.MutableData(
    name='delay_vars_data',
    value=data_df[delay_vars],
    dims=('date', 'delay_vars')
)

control_vars_data_ = pm.MutableData(
    name='control_vars_data',
    value=data_df[control_vars],
    dims=('date', 'control_vars')
)

fourier_vars_data_ = pm.MutableData(
    name='fourier_vars_data',
    value=data_df[fourier_vars],
    dims=('date', 'fourier_vars')
)

target_ = pm.MutableData(
    name='target',
    value=data_df[target_var],
    dims=('date')
)

# ------------------ Priors ------------------
# intercept = pm.HalfNormal('intercept', sigma=2)
sigma = pm.HalfNormal('sigma', sigma=2)

#### Delay vars
beta_delay_vars = pm.HalfNormal(
    name='beta_delay_vars',
    sigma=1,
    dims='delay_vars'
)

ls_delay_vars = pm.HalfNormal(
    name='ls_delay_vars',
    sigma=1,
    dims='delay_vars'
)

gp_delay_vars = pm.gp.Latent(
    cov_func=pm.gp.cov.ExpQuad(
        input_dim=len(delay_vars),
        # input_dim=1,
        ls=ls_delay_vars,
        # active_dims=[0, 1]
    ),
)

gp_prior_delay_vars = gp_delay_vars.prior(
    name='gp_prior_delay_vars',
    X=delay_vars_data_,
    # reparameterize=False,
    # dims=('date', 'delay_vars')
    dims=('delay_vars')
)

gamma_control_vars = pm.Normal(
    name='gamma_control_vars',
    mu=0,
    sigma=1,
    dims='control_vars'
)

ls_control_vars = pm.HalfNormal(
    name='ls_control_vars',
    sigma=1,
    dims='control_vars'
)

gp_control_vars = pm.gp.Latent(
    cov_func=pm.gp.cov.ExpQuad(
        input_dim=len(control_vars),
        # input_dim=1,
        ls=ls_control_vars,
        # active_dims=[2, 3]
    ),
)

gp_prior_control_vars = gp_control_vars.prior(
    name='gp_prior_control_vars',
    X=control_vars_data_,
    # reparameterize=False,
    dims=('control_vars')
)

phi_fourier_vars = pm.Laplace(
    name='phi_fourier_vars',
    mu=0,
    b=2,
    dims='fourier_vars'
)

ls_fourier_vars = pm.HalfNormal(
    name='ls_fourier_vars',
    sigma=1,
    dims='fourier_vars'
)

gp_fourier_vars = pm.gp.Latent(
    cov_func=pm.gp.cov.ExpQuad(
        input_dim=len(fourier_vars), 
        # input_dim=1,
        ls=ls_fourier_vars,
        # active_dims=[4, 5, 6, 7]
    ),
)

print(f"gp_fourier_vars.shape: {gp_fourier_vars})")

gp_prior_fourier_vars = gp_fourier_vars.prior(
    name='gp_prior_fourier_vars',
    X=fourier_vars_data_,
    # reparameterize=False,
    dims=('fourier_vars')
)


print(f"gp_prior_delay_vars.shape: {gp_prior_delay_vars.shape.eval()}")
print(f"gp_prior_control_vars.shape: {gp_prior_control_vars.shape.eval()}")
print(f"gp_prior_fourier_vars.shape: {gp_prior_fourier_vars.shape.eval()}")

delay_var_contributions = pm.Deterministic(
    name='delay_var_contribution',
    # var=pm.math.dot(gp_prior_delay_vars, delay_vars_data_),
    var=delay_vars_data_.T * gp_prior_delay_vars,
    # dims=('date', 'delay_vars')
    # dims=('delay_vars', 'date')
)
print(f"delay_var_contributions.shape: {delay_var_contributions.shape.eval()}")

control_var_contributions = pm.Deterministic(
    name='control_var_contribution',
    # var=pm.math.dot(gp_prior_control_vars, control_vars_data_),
    var=gp_prior_control_vars * control_vars_data_.T,
    # dims=('date', 'control_vars')
    dims=('control_vars', 'date')
)
print(f"control_var_contributions.shape: {control_var_contributions.shape.eval()}")

fourier_var_contributions = pm.Deterministic(
    name='fourier_var_contribution',
    # var=pm.math.dot(gp_prior_fourier_vars, fourier_vars_data_),
    var=gp_prior_fourier_vars * fourier_vars_data_.T,
    # dims=('date', 'fourier_vars')
    dims=('fourier_vars', 'date')
)
print(f"fourier_var_contributions.shape: {fourier_var_contributions.shape.eval()}")

print(delay_var_contributions.sum(axis=-1).shape.eval())

gp = delay_var_contributions.T.sum(axis=-1) +\
        control_var_contributions.T.sum(axis=-1) +\
            fourier_var_contributions.T.sum(axis=-1)



print(f"gp.shape: {gp.shape.eval()}")


likelihood = pm.Normal(
    name='likelihood',
    # mu=mu,
    mu=gp,
    sigma=sigma,
    observed=target_,
    dims='date'
)`  

Note: I'm using the same data as generated in the example: https://www.pymc-marketing.io/en/stable/notebooks/mmm/mmm_example.html

Note2: This approach doesn't even count with saturation and adstock transformations. This level of complexity was too much for me. I couldn't get to work even simpler versions...

I've tried multiple dims setting and multiple approaches which I found here:

• https://www.pymc-labs.io/blog-posts/spatial-gaussian-process-01/
• in your post: https://juanitorduz.github.io/gp_ts_pymc3/
• https://github.com/aloctavodia/BAP/blob/master/code/Chp7/07_Gaussian%20process.ipynb