A Fair Price for the Titanic? A Bayesian Analysis of the Price of Large Lego Sets

%matplotlib inline

import datetime
from functools import reduce

from aesara import shared, tensor as at
import arviz as az
from matplotlib import MatplotlibDeprecationWarning, pyplot as plt
from matplotlib.ticker import StrMethodFormatter
import numpy as np
import pandas as pd
import pymc3 as pm
from sklearn.preprocessing import StandardScaler
import seaborn as sns
import xarray as xr

FIGSIZE = (8, 6)

plt.rcParams['figure.figsize'] = FIGSIZE
sns.set(color_codes=True)

dollar_formatter = StrMethodFormatter("${x:,.0f}")

DATA_URL = 'https://austinrochford.com/resources/lego/brickset_19800101_20211098.csv.gz'

def to_datetime(year):
    return np.datetime64(f"{round(year)}-01-01")

COLS = [
    "Year released", "Set number",
    "Name", "Set type", "Theme", "Subtheme",
    "Pieces", "RRP$"
]

full_df = (pd.read_csv(DATA_URL, usecols=COLS)
             .dropna(subset=set(COLS) - {"Subtheme"}))
full_df["Year released"] = full_df["Year released"].apply(to_datetime)
full_df = (full_df.sort_values(["Year released", "Set number"])
                  .set_index("Set number"))

full_df.head()

full_df.tail()

full_df.describe()

CPI_URL = 'https://austinrochford.com/resources/lego/CPIAUCNS202100401.csv'

years = pd.date_range('1979-01-01', '2021-01-01', freq='Y') \
            + datetime.timedelta(days=1)
cpi_df = (pd.read_csv(CPI_URL, index_col="DATE", parse_dates=["DATE"])
            .loc[years])
cpi_df["to2021"] = cpi_df.loc["2021-01-01"] / cpi_df

full_df["RRP2021"] = (pd.merge(full_df, cpi_df,
                               left_on=["Year released"],
                               right_index=True)
                        .apply(lambda df: df["RRP$"] * df["to2021"],
                               axis=1))

FILTERS = [
    full_df["Set type"] == "Normal",
    full_df["Pieces"] > 10,
    full_df["Theme"] != "Duplo",
    full_df["Theme"] != "Service Packs",
    full_df["Theme"] != "Bulk Bricks",
    full_df["RRP2021"] > 0
]

df = full_df[reduce(np.logical_and, FILTERS)].copy()

df.head()

df.tail()

df.describe()

pieces = df["Pieces"].values
log_pieces = np.log(df["Pieces"].values)

scaler = StandardScaler().fit(log_pieces[:, np.newaxis])

def scale_log_pieces(log_pieces, scaler=scaler):
    return scaler.transform(log_pieces[:, np.newaxis])[:, 0]

std_log_pieces = scale_log_pieces(log_pieces)

rrp2021 = df["RRP2021"].values
log_rrp2021 = np.log(rrp2021)

theme_id, theme_map = df["Theme"].factorize(sort=True)

theme_mean_std_log_pieces = (pd.Series(std_log_pieces, index=df.index)
                               .groupby(df["Theme"])
                               .mean())

sub_id, sub_map = df["Subtheme"].factorize(sort=True)
n_sub = sub_map.size

sub_isnull = sub_id == -1

sub_mean_std_log_pieces = (pd.Series(std_log_pieces, index=df.index)
                             .groupby(df["Subtheme"])
                             .mean())

year = (df["Year released"]
          .dt.year
          .values)
t = year - year.min()

def hierarchical_normal(name, *, dims, μ=None):
    if μ is None:
        μ = pm.Normal(f"μ_{name}", 0., 2.5)

    Δ = pm.Normal(f"Δ_{name}", 0., 1., dims=dims)
    σ = pm.HalfNormal(f"σ_{name}", 2.5)
    
    return pm.Deterministic(name, μ + Δ * σ, dims=dims)

def hierarchical_normal_with_mean(name, x_mean, *, dims,
                                  μ=None, mean_name="γ"):
    mean_coef = pm.Normal(f"{mean_name}_{name}", 0., 2.5)
        
    return pm.Deterministic(
        name,
        hierarchical_normal(f"_{name}", dims=dims, μ=μ) \
            + mean_coef * x_mean,
        dims=dims
    )

def gaussian_random_walk(name, *, dims, innov_scale=1.):
    Δ = pm.Normal(f"Δ_{name}", 0., innov_scale,  dims=dims)

    return pm.Deterministic(name, Δ.cumsum(), dims=dims)

coords = {
    "set": df.index.values,
    "sub": sub_map,
    "theme": theme_map,
    "year": np.unique(year)
}

with pm.Model(coords=coords) as model:
    β0_0 = pm.Normal("β0_0", 0., 2.5)
    β0_t = gaussian_random_walk("β0_t", dims="year",
                                innov_scale=0.1)
    
    β0_theme = hierarchical_normal_with_mean(
        "β0_theme",
        theme_mean_std_log_pieces.values,
        dims="theme", μ=0.
    )
    
    β0_sub_null = pm.Normal("β0_sub_null", 0., 2.5)
    β0_sub_nn = hierarchical_normal_with_mean(
        "β0_sub_nn",
        sub_mean_std_log_pieces.values,
        dims="sub", μ=0., mean_name="λ"
    )
    β0_sub_i = at.switch(
        sub_isnull,
        β0_sub_null,
        β0_sub_nn[at.clip(sub_id, 0, n_sub)]
    )
    
    β0_i = β0_0 + β0_t[t] + β0_theme[theme_id] + β0_sub_i

with model:
    β_pieces_0 = pm.Normal("β_pieces_0", 0., 2.5)
    
    β_pieces_theme = hierarchical_normal_with_mean(
        "β_pieces_theme",
        theme_mean_std_log_pieces.values,
        dims="theme", μ=0.
    )
    
    β_pieces_sub_null = pm.Normal("β_pieces_sub_null", 0., 2.5)
    β_pieces_sub_nn = hierarchical_normal_with_mean(
        "β_pieces_sub_nn",
        sub_mean_std_log_pieces.values,
        dims="sub", μ=0., mean_name="λ"
    )
    β_pieces_sub_i = at.switch(
        sub_isnull,
        β_pieces_sub_null,
        β_pieces_sub_nn[at.clip(sub_id, 0, n_sub)]
    )

    β_pieces_i = β_pieces_0 + β_pieces_theme[theme_id] + β_pieces_sub_i

with model:
    σ = pm.HalfNormal("σ", 5.)
    μ = β0_i + β_pieces_i * std_log_pieces - 0.5 * σ**2
    
    obs = pm.Normal("obs", μ, σ, dims="set",
                    observed=log_rrp2021)

CHAINS = 3
SEED = 123456789

SAMPLE_KWARGS = {
    'cores': CHAINS,
    'random_seed': [SEED + i for i in range(CHAINS)],
    'return_inferencedata': True
}

with model:
    trace = pm.sample(**SAMPLE_KWARGS)
    pp_trace = pm.to_inference_data(
        posterior_predictive=pm.sample_posterior_predictive(trace)
    )

trace.extend(pp_trace)