Hello, PyMC community!

Please help with the following question.

I have an e-commerce store working for two cities (call them city A and city B). Both cities share a very common pattern of people buying things, so the two timelines resemble each other in terms of trend and fluctuations (please have a look at the screenshot attached).

I believe that both cities share some **common latent dynamic of purchasing behaviour** and would like to be able to extract it.

My attempt to model it is to fit the two cities within the model.

- both cities share the same trend
- both cities share the same weekly seasonality*
- but have different intercepts

*(the pattern is the same, but I can see that itâs proportionally a bit different, so I add a âsqueezerâ, to account for proportional changes).

Hereâs the specification

```
day_of_week_idx, day_of_week = df["day_of_week"].factorize(sort=True)
date = df["date_day"].to_numpy()
coords = {
"day_of_week": day_of_week,
"date" :date
}
with pm.Model(coords=coords) as model:
# target variables
sales_city_a = pm.MutableData(name="sales_city_a", value=df['sales_city_a'])
sales_city_b = pm.MutableData(name="sales_city_b", value=df['sales_city_b'])
# Intercepts: different for both cities
intercept_city_a = pm.Normal(name="intercept_city_a", mu=0, sigma=200)
intercept_city_b = pm.Normal(name="intercept_city_b", mu=0, sigma=200)
# trend: common for both cities
general_trend = pm.Normal(name="general_trend", mu=150, sigma=200)
trend = pm.MutableData(name="trend", value=df['trend'])
# Weekly seasonality - common for both cities, but in one city is proportional to the other, so we use 'squeeze coefficient'
b_day_of_week = pm.ZeroSumNormal(name="b_day_of_week", sigma=100, dims="day_of_week")
day_of_week = pm.MutableData(name="day_of_week", value=df['day_of_week'])
beta_week_squeeze = pm.Normal(name="beta_week_squeeze", mu=0, sigma=2)
# Modeled behaviour of two cities
mu_city_a = pm.Deterministic(name="mu_city_a", var=intercept_city_a + trend*general_trend + b_day_of_week[day_of_week_idx]*beta_week_squeeze*day_of_week)
mu_city_b = pm.Deterministic(name="mu_city_b", var=intercept_city_b + trend*general_trend + b_day_of_week[day_of_week_idx]*day_of_week)
sigma = pm.HalfNormal(name="sigma", sigma=200)
# fit two cities
pm.Normal(name="l_city_a", mu=mu_city_a, sigma=sigma, observed=sales_city_a)
pm.Normal(name="l_city_b", mu=mu_city_b, sigma=sigma, observed=ftr_ios)
```

As a result I get a model that fits the data quite well.

However, there are some days when something happens in both cities. So the purchasing pattern is still shared across the two cities, but is not captured by the model (see screenshot below).

**So my question is:**

What is the best way to tell the model about those cases? I would like to **capture not all deviations, but only those that are shared by both cities simultaneously**. One option is to mark all those days as dummy variables. Or is it correct to collect all the âerrorsâ to some separate variable? Or can those be modelled as some process (like Gaussian Random Walk, for example?)

**Why am I doing that?**

We have been running multiple **advertising activities** in the two cities at different times. So I assume that each cityâs timeline can be decomposed to the purchasing pattern (common for both cities) + the impact of advertising.

Therefore, Iâd like to extract the common purchasing pattern as accurately as possible, so that the model does not confuse it with the impact of advertising.

Appreciate any recommendations!

**PS.** When you want to model the impact of advertising on your sales, one simpler thing you can do is just put sales in a city as target variable and put advertising activities as predictors. However, in this case the model over/underestimates the impact, often confusing it with other external factors (which are out of our control and often immeasurable). So we want to take advantage of both cities sharing the same pattern and assume that this common pattern would capture all the external factors.