Bayesian piecewise ITS of arXiv Computer Science submission growth

arxiv_cs_piecewise_its.py

"""Bayesian segmented analysis of arXiv CS primary submissions."""

import re
from urllib.request import urlopen

import causalpy as cp
import matplotlib.pyplot as plt
import pandas as pd

MONTHLY_COUNTS = re.compile(
    r"<b>(\d+)</b> \+ \d+ \([A-Z][a-z]{2} \d{4}\)"
)


def primary_cs_submissions() -> pd.DataFrame:
    """Download primary (not cross-listed) CS counts from arXiv's yearly pages."""
    rows = []
    for year in range(2015, pd.Timestamp.today().year + 1):
        html = urlopen(f"https://arxiv.org/year/cs/{year}", timeout=30).read().decode()
        rows.extend(
            (f"{year}-{month:02d}-01", int(count))
            for month, count in enumerate(MONTHLY_COUNTS.findall(html), start=1)
        )
    data = pd.DataFrame(rows, columns=["date", "submissions"])
    data["date"] = pd.to_datetime(data["date"])
    data = data[
        data.date < pd.Timestamp.today().replace(day=1).normalize()
    ].copy()
    assert pd.DatetimeIndex(data.date).equals(
        pd.date_range(data.date.min(), data.date.max(), freq="MS")
    )
    data["t"] = (data.date - data.date.min()).dt.days / 365.25
    data["y"] = data.submissions / 1_000
    return data


data = primary_cs_submissions()
result = cp.PiecewiseITS(
    data,
    formula=(
        "y ~ 1 + t + I(ramp(date, '2020-05-01') / 365.25)"
        " + I(ramp(date, '2023-01-01') / 365.25)"
        " + I(ramp(date, '2026-01-01') / 365.25)"
    ),
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"random_seed": 42, "progressbar": False}
    ),
)
fig, axes = result.plot(show=False)
axes[0].set_ylabel("submissions (thousands)")
fig.savefig("arxiv-cs-piecewise-its.png", dpi=200, bbox_inches="tight")

beta = result.idata.posterior["beta"].sel(treated_units="unit_0")
slopes = [beta.sel(coeffs="t")]
for ramp in (label for label in result.labels if "ramp(" in label):
    slopes.append(slopes[-1] + beta.sel(coeffs=ramp))

slope_fig, slope_ax = plt.subplots(figsize=(8, 4))
slope_ax.violinplot(
    [slope.values.ravel() for slope in slopes], showmeans=True, showextrema=False
)
slope_ax.axhline(0, color="black", linestyle="--", linewidth=1)
slope_ax.set(
    xticks=range(1, 5),
    xticklabels=["2015–Apr 2020", "May 2020–Dec 2022", "Jan 2023–Dec 2025", "Jan 2026–"],
    ylabel="slope (thousand submissions/year)",
    title="Posterior trend slopes",
)
slope_fig.tight_layout()
slope_fig.savefig("arxiv-cs-piecewise-its-slopes.png", dpi=200, bbox_inches="tight")