| 769 | \pending{Rao--Stirling and Leinster--Cobbold indices per category and period, not yet computed} |
769 | \pending{Rao--Stirling and Leinster--Cobbold indices per category and period, not yet computed} |
| 770 | until the computation is added and audited. |
770 | until the computation is added and audited. |
| 771 | |
771 | |
| |
772 | +\section{Citation and collaboration networks} |
| |
773 | +\label{sec:m-networks} |
| |
774 | + |
| |
775 | +\subsection{The citation graph and its boundary} |
| |
776 | +\label{sec:m-citgraph} |
| |
777 | + |
| |
778 | +Citation edges come from the OpenAlex \texttt{referenced\_works} lists, which |
| |
779 | +are verified by construction in the sense that each edge is an assertion of the |
| |
780 | +source database rather than an inference of this study. The script |
| |
781 | +\texttt{knowledge\_graph/01\_citation\_graph.py} splits the references of |
| |
782 | +eligible works into two edge sets: |
| |
783 | +\begin{align} |
| |
784 | + \mathcal{E}_{\mathrm{int}} &= \{ (w, v) : w \in \mathcal{W},\ v \in \mathcal{W},\ v \in \mathrm{refs}(w) \},\\ |
| |
785 | + \mathcal{E}_{\mathrm{ext}} &= \{ (w, v) : w \in \mathcal{W},\ v \notin \mathcal{W},\ v \in \mathrm{refs}(w),\ v \in \text{slim table} \}, |
| |
786 | +\end{align} |
| |
787 | +the internal edges forming the citation graph of computer science proper and |
| |
788 | +the external edges recording what the corpus cites outside itself, with the |
| |
789 | +cited work's publication year, type and primary OpenAlex field resolved from the |
| |
790 | +slim table of all works. References that resolve to no OpenAlex record are |
| |
791 | +counted (the declared reference count of each work is stored) but cannot be |
| |
792 | +typed. The yearly tables derived from these sets are the inflows by cited |
| |
793 | +field (how much of each citing year's references go to mathematics, medicine, |
| |
794 | +the social sciences, and so on), the distribution of citation age |
| |
795 | +$y(w) - y(v)$ by citing year, and the reference coverage (declared, internal, |
| |
796 | +external) by year. Coverage matters for every network statistic that follows: |
| |
797 | +if the share of resolvable references rises over time, as it does in OpenAlex |
| |
798 | +because recent works are better matched, then internal edge counts rise for |
| |
799 | +reasons unrelated to citation behaviour, and all flow statistics are therefore |
| |
800 | +expressed as shares of resolved edges within the citing year. |
| |
801 | + |
| |
802 | +\subsection{Category flow matrices} |
| |
803 | + |
| |
804 | +For each citing year $t$ the script \texttt{knowledge\_graph/02\_category\_flows.py} |
| |
805 | +forms the $K \times K$ matrix |
| |
806 | +\begin{equation} |
| |
807 | + F_t[i, j] \;=\; \bigl|\{ (w, v) \in \mathcal{E}_{\mathrm{int}} : y(w) = t,\ \kappa^{\ast}(w) = i,\ \kappa^{\ast}(v) = j \}\bigr| , |
| |
808 | + \label{eq:flow} |
| |
809 | +\end{equation} |
| |
810 | +the number of citations made in year $t$ by works whose primary category is $i$ |
| |
811 | +to works whose primary category is $j$, together with the row-normalised |
| |
812 | +shares $F_t[i,j] / \sum_{j'} F_t[i,j']$. Works without a primary category |
| |
813 | +contribute no edges to $F_t$; their references are retained in |
| |
814 | +$\mathcal{E}_{\mathrm{int}}$ but fall outside the matrix. Aggregating rows and |
| |
815 | +columns by area through $a(\cdot)$ gives the $8 \times 8$ area flow matrix, and |
| |
816 | +the cross-area share of citing area $\alpha$ in year $t$ is |
| |
817 | +\begin{equation} |
| |
818 | + X_{\alpha, t} \;=\; \frac{\sum_{i : a(i)=\alpha}\ \sum_{j : a(j) \ne \alpha} F_t[i,j]}{\sum_{i : a(i)=\alpha}\ \sum_{j} F_t[i,j]}, |
| |
819 | + \label{eq:cross-area} |
| |
820 | +\end{equation} |
| |
821 | +the fraction of an area's internal references that leave the area. This is the |
| |
822 | +statistic behind hypothesis H4. It is not a measure of interdisciplinarity in |
| |
823 | +the sense of Section~\ref{sec:m-diversity}, because it treats all other areas |
| |
824 | +as equally distant; it is a boundary-crossing rate, and its interpretation |
| |
825 | +depends on the stability of the boundaries, which the classifier fixes once |
| |
826 | +for all years. A co-classification matrix counts, for each year, the works |
| |
827 | +carrying both category $i$ and category $j$ in their full label set $C(w)$; it |
| |
828 | +is reported alongside the flows but, as noted in |
| |
829 | +Section~\ref{sec:m-classification}, it is sensitive to abstract length and |
| |
830 | +vocabulary in a way the primary-label flows are not. |
| |
831 | + |
| |
832 | +\subsection{Modularity of the yearly flow network} |
| |
833 | + |
| |
834 | +The degree to which citations stay within categories, or within areas, is |
| |
835 | +summarised by the modularity of the weighted directed graph with adjacency |
| |
836 | +$W = F_t$, evaluated at the fixed partition given by the categories themselves |
| |
837 | +(each category its own community) and at the coarser partition given by the |
| |
838 | +areas. The implemented statistic is |
| |
839 | +\begin{equation} |
| |
840 | + Q(W, g) \;=\; \frac{1}{m} \sum_{i,j} \Bigl( W_{ij} - \frac{k^{\mathrm{out}}_i\, k^{\mathrm{in}}_j}{m} \Bigr)\, \mathbf{1}[g_i = g_j], |
| |
841 | + \qquad m = \sum_{i,j} W_{ij},\quad k^{\mathrm{out}}_i = \sum_j W_{ij},\quad k^{\mathrm{in}}_j = \sum_i W_{ij}, |
| |
842 | + \label{eq:modularity} |
| |
843 | +\end{equation} |
| |
844 | +which is the directed generalisation by \citet{leicht2008directed} of the |
| |
845 | +Newman--Girvan modularity \citep{newman2004community}: the observed weight |
| |
846 | +within communities minus the weight expected if each citation were rewired at |
| |
847 | +random while preserving every node's out-strength and in-strength. Because the |
| |
848 | +partition is given rather than optimised, $Q$ is used here as a descriptive |
| |
849 | +index of assortativity by category, not as the objective of a community |
| |
850 | +detection algorithm; the resolution limit that affects optimised modularity |
| |
851 | +\citep{fortunato2007resolution, fortunato2010community} therefore does not |
| |
852 | +arise, but the index inherits the dependence of the flow matrix on the |
| |
853 | +classifier. Two values are reported per year, $Q_{\mathrm{cat}}$ with |
| |
854 | +$g_i = i$ and $Q_{\mathrm{area}}$ with $g_i = a(i)$, together with the self-share |
| |
855 | +$\operatorname{tr}(F_t)/m$, which is the uncorrected fraction of citations |
| |
856 | +staying within the citing category. Hypothesis H4 predicts a decrease in |
| |
857 | +$Q_{\mathrm{area}}$ in the post period. A fall in $Q$ can also arise if the |
| |
858 | +composition of citing works shifts toward categories that were always |
| |
859 | +boundary-crossing (large language models citing natural-language processing |
| |
860 | +and machine-learning theory), so Part IX reports $Q$ together with the |
| |
861 | +composition-adjusted cross-area shares obtained by fixing the row weights at |
| |
862 | +their baseline values. |
| |
863 | + |
| |
864 | +\subsection{Co-authorship, institution and country networks} |
| |
865 | + |
| |
866 | +The plan lists author, institution and country collaboration networks built |
| |
867 | +with \texttt{igraph}. The implemented collaboration layer is tabular rather than |
| |
868 | +graph-analytic: country pairs with at least fifty joint works per period |
| |
869 | +(edges of the country network, with their weights), institution counts by |
| |
870 | +period under full and fractional counting, and institution-type shares. The |
| |
871 | +country network is thus available as an edge list from which the figures of |
| |
872 | +Part VIII are drawn, and degree, strength and density can be computed from it; |
| |
873 | +author-level co-authorship networks are not constructed in the current |
| |
874 | +pipeline, because author disambiguation in OpenAlex is of uneven quality over |
| |
875 | +the period and an author graph of this size would be dominated by that noise. |
| |
876 | +This is a narrowing of the plan, recorded here; the country and institution |
| |
877 | +tables are the collaboration estimands actually reported. |
| |
878 | + |
| |
879 | +\section{Semantic mapping} |
| |
880 | +\label{sec:m-semantic} |
| |
881 | + |
| |
882 | +\subsection{Document embeddings} |
| |
883 | + |
| |
884 | +Each eligible work with an abstract is represented by a dense vector computed |
| |
885 | +by a scientific sentence-embedding model from the string |
| |
886 | +``$\text{title}$. $\text{abstract}[1{:}2000]$'' (the abstract truncated to |
| |
887 | +2,000 characters), tokenised to at most 256 word pieces, and L2-normalised, so |
| |
888 | +that the inner product of two vectors is their cosine similarity |
| |
889 | +(\texttt{embeddings/embed\_corpus.py}). The primary model is SPECTER2 base |
| |
890 | +\citep{singh2023scirepeval}, the successor of SPECTER |
| |
891 | +\citep{cohan2020specter}, a transformer trained with a citation-informed |
| |
892 | +triplet objective so that papers that cite one another are close in the |
| |
893 | +embedding space; this inductive bias is the appropriate one for mapping a |
| |
894 | +literature by its intellectual proximity rather than by lexical overlap. The |
| |
895 | +protocol requires a second embedding family for comparison; the benchmark |
| |
896 | +script evaluates SciNCL \citep{ostendorff2022scincl}, a general-purpose |
| |
897 | +BGE model and a MiniLM sentence-transformer \citep{reimers2019sbert} for |
| |
898 | +throughput and dimensionality on the same hardware, and the atlas script |
| |
899 | +accepts any of them through its model argument, so the comparison is a |
| |
900 | +re-run with a different tag. The embeddings are stored as float16 shards with |
| |
901 | +their work identifiers and a per-host manifest (model name, maximum sequence |
| |
902 | +length, device, timings), computed on several nodes of the cluster in |
| |
903 | +parallel; float16 storage changes cosine similarities by less than |
| |
904 | +$10^{-3}$, which is below any threshold used downstream. |
| |
905 | + |
| |
906 | +Three assumptions underlie everything built on these vectors. The abstract is |
| |
907 | +assumed to be a faithful summary of the work's content, which it is for most |
| |
908 | +research articles and is not for extended abstracts and some proceedings |
| |
909 | +papers. The embedding geometry is assumed stable across the period: a model |
| |
910 | +trained on literature up to its release date represents 2026 work in a space |
| |
911 | +shaped by earlier work, so genuinely new vocabulary is embedded by analogy to |
| |
912 | +what the model knows, and the novelty statistic below measures distance in |
| |
913 | +that space, not in some ideal space of ideas. And cosine similarity is assumed |
| |
914 | +to be a sufficient proxy for topical proximity; this is the standard |
| |
915 | +assumption of the field and is supported by the retrieval benchmarks on which |
| |
916 | +these models are evaluated, with the caveat that all such benchmarks are |
| |
917 | +dominated by pre-2023 literature. |
| |
918 | + |
| |
919 | +\subsection{Atlas: projection and clustering} |
| |
920 | + |
| |
921 | +The two-dimensional atlas and the cluster structure are computed on a uniform |
| |
922 | +random sample of 200,000 eligible works with embeddings (seed 44), because |
| |
923 | +density-based clustering of the full corpus in the native space is neither |
| |
924 | +computationally practical nor, given the noise level of abstracts, more |
| |
925 | +informative. UMAP \citep{mcinnes2018umap} with the cosine metric, 30 nearest |
| |
926 | +neighbours and a minimum distance of zero projects the sample to 50 |
| |
927 | +dimensions; HDBSCAN \citep{campello2013hdbscan, mcinnes2017hdbscan} with a |
| |
928 | +minimum cluster size of 200 and a minimum sample count of 20 is run on the |
| |
929 | +50-dimensional projection; a second UMAP to two dimensions (minimum distance |
| |
930 | +0.1) is used for display only. The 50-dimensional intermediate step is what |
| |
931 | +makes density-based clustering workable on 768-dimensional normalised vectors, |
| |
932 | +whose pairwise distances concentrate; it is also a source of distortion, since |
| |
933 | +UMAP preserves local neighbourhoods and not global distances, so cluster |
| |
934 | +\emph{membership} is meaningful and the distances between clusters on the |
| |
935 | +two-dimensional map are not. HDBSCAN labels points in low-density regions as |
| |
936 | +noise ($-1$), and the noise share is reported with the cluster count. Each |
| |
937 | +cluster is described by the twelve title terms with the highest mean TF--IDF |
| |
938 | +weight (unigrams and bigrams, English stop words removed, minimum document |
| |
939 | +frequency five), by its modal primary category and that category's share, |
| |
940 | +and by the share of its members in the post period; these descriptors are |
| |
941 | +how clusters are related to the ontology, and the comparison of the two |
| |
942 | +partitions (bottom-up clusters against top-down categories) is one of the |
| |
943 | +checks on the ontology reported in Part IX. |
| |
944 | + |
| |
945 | +\subsection{Stability} |
| |
946 | + |
| |
947 | +Cluster stability is measured by a split-half design. Five times, the sample |
| |
948 | +is randomly permuted and divided into two halves; HDBSCAN (minimum cluster |
| |
949 | +size 100, since each half has half the points) is run separately on the two |
| |
950 | +halves of the \emph{same} 50-dimensional projection; a 5-nearest-neighbour |
| |
951 | +classifier trained on the first half's labels predicts labels for the second |
| |
952 | +half; and the adjusted Rand index \citep{hubert1985comparing} |
| |
953 | +\begin{equation} |
| |
954 | + \mathrm{ARI} \;=\; \frac{\sum_{ij}\binom{n_{ij}}{2} - \bigl[\sum_i \binom{a_i}{2}\sum_j \binom{b_j}{2}\bigr]/\binom{n}{2}} |
| |
955 | + {\tfrac12\bigl[\sum_i \binom{a_i}{2} + \sum_j \binom{b_j}{2}\bigr] - \bigl[\sum_i \binom{a_i}{2}\sum_j \binom{b_j}{2}\bigr]/\binom{n}{2}} |
| |
956 | +\end{equation} |
| |
957 | +between the transferred labels and the second half's own labels is recorded |
| |
958 | +(here $n_{ij}$ is the contingency table of the two labellings and $a_i$, $b_j$ |
| |
959 | +its margins). The mean of the five values is the reported stability. Two |
| |
960 | +caveats bound what this measures. The projection is fitted once on the full |
| |
961 | +sample, so the design tests the stability of the clustering conditional on the |
| |
962 | +projection and not the stability of UMAP itself, whose stochastic optimisation |
| |
963 | +is fixed by the seed; the protocol's phrase ``bootstrap ARI'' describes this |
| |
964 | +split-half procedure, which resamples without replacement. And the |
| |
965 | +nearest-neighbour transfer is itself a smoothing step: it maps noise points of |
| |
966 | +the first half to the majority label of their neighbours, which raises the |
| |
967 | +index relative to a direct comparison of two independent clusterings of the |
| |
968 | +same points. |
| |
969 | + |
| |
970 | +\subsection{Category centroids, drift and novelty} |
| |
971 | +\label{sec:m-novelty} |
| |
972 | + |
| |
973 | +In the native embedding space, the centroid of category $k$ in year $y$ is |
| |
974 | +the normalised mean vector |
| |
975 | +\begin{equation} |
| |
976 | + \mu_{k,y} \;=\; \frac{\sum_{w : \kappa^{\ast}(w)=k,\ y(w)=y} e_w}{\bigl\|\sum_{w : \kappa^{\ast}(w)=k,\ y(w)=y} e_w\bigr\|}, |
| |
977 | +\end{equation} |
| |
978 | +computed only for cells with at least 100 works, and the year-over-year drift |
| |
979 | +of a category is the cosine distance $1 - \mu_{k,y}^{\top}\mu_{k,y-1}$. Period |
| |
980 | +centroids average the yearly centroids over calendar years 2018--2022 |
| |
981 | +(baseline) and 2023--2026 (post), and the inter-category cosine similarity |
| |
982 | +matrices of the two periods are the semantic counterpart of hypotheses H4 and |
| |
983 | +H5: convergence between two categories appears as an increase in the cosine |
| |
984 | +similarity of their centroids. The period boundary used here is the calendar |
| |
985 | +year, not the 30~November 2022 date of Section~\ref{sec:m-notation}; the |
| |
986 | +eleven months of 2022 that belong to the baseline are kept in the baseline |
| |
987 | +and December 2022 is attached to it rather than to the post period, a |
| |
988 | +deliberate simplification at the level of yearly centroids that moves a |
| |
989 | +single month across the boundary. |
| |
990 | + |
| |
991 | +The novelty of a work is its distance to the recent past of its own category: |
| |
992 | +\begin{equation} |
| |
993 | + \nu(w) \;=\; 1 - e_w^{\top}\,\bar\mu_{\kappa^{\ast}(w),\ m(w)-24\,:\,m(w)-1}, \qquad |
| |
994 | + \bar\mu_{k,\,a:b} = \frac{\sum_{t=a}^{b} \sum_{v : \kappa^{\ast}(v)=k,\ m(v)=t} e_v}{\bigl\|\,\cdot\,\bigr\|}, |
| |
995 | + \label{eq:novelty} |
| |
996 | +\end{equation} |
| |
997 | +the cosine distance between the work's vector and the normalised sum of the |
| |
998 | +vectors of all works of the same primary category whose first month lies in |
| |
999 | +the 24 months strictly before the work's own month. A work in the first months |
| |
1000 | +of a category, with no trailing window, has no novelty value. The |
| |
1001 | +definition uses only information dated before the work, which is what allows |
| |
1002 | +it to serve as a signal in the out-of-time detector of |
| |
1003 | +Section~\ref{sec:m-emerging}; but the centroid is computed with the |
| |
1004 | +\emph{final} primary categories of the earlier works, assigned by rules fixed |
| |
1005 | +in 2026, so the category boundaries themselves are not out-of-time. This is a |
| |
1006 | +limitation shared by every category-conditional statistic in the paper and |
| |
1007 | +is discussed in Section~\ref{sec:m-emerging}. Novelty defined as distance to a |
| |
1008 | +category centroid rewards works that sit between categories as much as works |
| |
1009 | +that are new; the interpretation of high novelty as ``new'' rather than |
| |
1010 | +``interdisciplinary'' requires the clusters and the flow statistics as |
| |
1011 | +corroboration, and the literature on novelty indicators warns that |
| |
1012 | +citation-based and text-based novelty measures do not agree well |
| |
1013 | +\citep{uzzi2013atypical, wang2017novelty}. |
| |
1014 | + |
| |
1015 | +\section{Structural change in monthly series} |
| |
1016 | +\label{sec:m-structural} |
| |
1017 | + |
| |
1018 | +The question behind RQ2 and hypotheses H1 and H3 is whether the monthly |
| |
1019 | +series of computer-science output, in total and by category, changed level or |
| |
1020 | +slope around the reference date, and whether any such change stands out |
| |
1021 | +against the ordinary variability of the series. Four methods are applied to |
| |
1022 | +every series by \texttt{econometrics/run\_structural\_change.py}, using the |
| |
1023 | +pure functions of \texttt{econometrics/structural\_change.py} (unit-tested on a |
| |
1024 | +simulated series with a known break). None of them identifies a cause. An |
| |
1025 | +interrupted time series with a single pre-specified date answers ``did the |
| |
1026 | +series change at that date, net of trend and season?''; the placebo |
| |
1027 | +distribution answers ``is that date special?''; the multiple-break and |
| |
1028 | +changepoint methods answer ``where does the series change, if one does not |
| |
1029 | +specify a date?''; and the Bayesian online detector answers ``when would an |
| |
1030 | +observer following the series month by month have concluded that it had |
| |
1031 | +changed?''. The reference date is a historical marker, and the manuscript |
| |
1032 | +uses no causal vocabulary about it (CLAUDE.md §6.5). |
| |
1033 | + |
| |
1034 | +\subsection{Series construction and truncation} |
| |
1035 | +\label{sec:m-truncation} |
| |
1036 | + |
| |
1037 | +Each series is a vector of monthly counts $N_t$ on the month-precision |
| |
1038 | +subset, from February 2018 to the last complete month $t_{\mathrm{L}}$, by |
| |
1039 | +default June 2026, with sensitivity runs at March and August 2026. The first |
| |
1040 | +month of the window is dropped because January 2018 is contaminated by the |
| |
1041 | +year-precision convention of Section~\ref{sec:m-notation} in the other |
| |
1042 | +direction (year-precision works with a genuine January date are not |
| |
1043 | +recoverable). Category series $N_{k,t}$ use the multi-label count restricted to |
| |
1044 | +month-precision works, and a series is analysed only if it has at least 48 |
| |
1045 | +non-zero months and at least 500 works in total; smaller series yield |
| |
1046 | +unstable fits and are reported descriptively. All methods operate on |
| |
1047 | +$y_t = \log(N_t + 1)$, the shift making the transform defined for empty months |
| |
1048 | +and immaterial for counts in the hundreds or thousands. |
| |
1049 | + |
| |
1050 | +Right-truncation is handled by truncation, not by correction. The plan |
| |
1051 | +proposed to correct recent months by the indexing-delay profile estimated from |
| |
1052 | +the creation dates; the implemented pipeline estimates that profile |
| |
1053 | +(\texttt{indexing\_lag.csv}: for each publication month, the median and 90th |
| |
1054 | +percentile of $d_{\mathrm{cr}}(w) - d_{\mathrm{pub}}(w)$ and the shares indexed |
| |
1055 | +within 90 and after 180 days) and reports it, but the structural-change runs |
| |
1056 | +simply stop the series at $t_{\mathrm{L}}$ and vary $t_{\mathrm{L}}$. The |
| |
1057 | +reason is that an inflation correction based on the lag profile of earlier |
| |
1058 | +years assumes that the profile is stationary, and the profile is itself |
| |
1059 | +affected by the events under study (OpenAlex's own ingestion changed over the |
| |
1060 | +period); a correction would import that assumption into every test. Truncating |
| |
1061 | +to the last month at which the lag profile indicates near-complete indexing |
| |
1062 | +loses a few months of data but introduces no model. Donut windows (dropping |
| |
1063 | +the months adjacent to the event) are not implemented; the alternative event |
| |
1064 | +dates of Section~\ref{sec:m-altdates} serve the related purpose of showing |
| |
1065 | +whether the estimated change is tied to a particular month. |
| |
1066 | + |
| |
1067 | +\subsection{Interrupted time series with seasonality and HAC errors} |
| |
1068 | +\label{sec:m-its} |
| |
1069 | + |
| |
1070 | +With $t = 0, 1, \dots, n-1$ indexing the months and $t_e$ the index of the |
| |
1071 | +event month, the segmented regression is |
| |
1072 | +\begin{equation} |
| |
1073 | + y_t \;=\; \beta_0 + \beta_1\, t + \beta_2\, \mathbf{1}[t \ge t_e] + \beta_3\, (t - t_e)\,\mathbf{1}[t \ge t_e] |
| |
1074 | + + \sum_{j=2}^{12} \gamma_j\, \mathbf{1}[\mathrm{month}(t) = j] + \varepsilon_t , |
| |
1075 | + \label{eq:its} |
| |
1076 | +\end{equation} |
| |
1077 | +the standard interrupted-time-series specification \citep{bernal2017its} with |
| |
1078 | +eleven month-of-year dummies (January as reference). The coefficient |
| |
1079 | +$\beta_2$ is the change in log level at the event, $\beta_3$ the change in |
| |
1080 | +monthly log slope, and $\beta_1$ the pre-event slope; since $y$ is a log count, |
| |
1081 | +$e^{\beta_2} - 1$ is the proportional level shift and $\beta_1$ and $\beta_3$ |
| |
1082 | +are monthly growth rates. The event month is coded as the first full month |
| |
1083 | +after the reference date, December 2022 for the 30~November 2022 marker, and |
| |
1084 | +likewise for the alternative dates. Ordinary least squares estimates the |
| |
1085 | +coefficients; the covariance matrix is the heteroskedasticity- and |
| |
1086 | +autocorrelation-consistent estimator of \citet{newey1987hac} with a Bartlett |
| |
1087 | +kernel and bandwidth six months, |
| |
1088 | +\begin{equation} |
| |
1089 | + \hat V = (X^{\top}X)^{-1} \Bigl( \hat\Gamma_0 + \sum_{\ell=1}^{6} \bigl(1 - \tfrac{\ell}{7}\bigr)(\hat\Gamma_\ell + \hat\Gamma_\ell^{\top}) \Bigr) (X^{\top}X)^{-1}, |
| |
1090 | + \qquad \hat\Gamma_\ell = \sum_{t} \hat\varepsilon_t \hat\varepsilon_{t-\ell}\, x_t x_{t-\ell}^{\top}, |
| |
1091 | +\end{equation} |
| |
1092 | +which makes the reported standard errors robust to serial correlation of the |
| |
1093 | +residuals up to half a year. The bandwidth is fixed rather than chosen by the |
| |
1094 | +data-dependent rule of \citet{andrews1991hac}; with roughly one hundred |
| |
1095 | +observations the difference is immaterial and a fixed choice is easier to |
| |
1096 | +preregister. The fitted values are back-transformed as $e^{\hat y_t} - 1$ for |
| |
1097 | +the figures. The specification's known weaknesses are those of any single-break |
| |
1098 | +segmented regression: it assumes the change is immediate and permanent, it |
| |
1099 | +attributes to $\beta_2$ any curvature of the pre-trend that a linear |
| |
1100 | +$\beta_1 t$ cannot absorb, and with a reference date placed near the centre of |
| |
1101 | +the sample it is well identified only because the post period is long. The |
| |
1102 | +placebo analysis is the main guard against the second weakness. |
| |
1103 | + |
| |
1104 | +\subsection{Placebo distribution} |
| |
1105 | +\label{sec:m-placebo} |
| |
1106 | + |
| |
1107 | +The ITS of equation~\eqref{eq:its} is re-estimated with the event placed at |
| |
1108 | +every month $t'$ from the 25th month of the series to the 12th month before |
| |
1109 | +its end, yielding the $t$-statistics $\{ t^{(2)}_{t'} \}$ and $\{ t^{(3)}_{t'} \}$ |
| |
1110 | +of the level and slope changes as a function of the hypothetical event date. |
| |
1111 | +The placebo $p$-value of the true event is its rank among these: |
| |
1112 | +\begin{equation} |
| |
1113 | + p^{\mathrm{plac}}_{\beta_2} \;=\; \frac{1}{|\mathcal{T}'|}\sum_{t' \in \mathcal{T}'} \mathbf{1}\bigl[\, |t^{(2)}_{t'}| \ge |t^{(2)}_{t_e}| \,\bigr], |
| |
1114 | + \label{eq:placebo} |
| |
1115 | +\end{equation} |
| |
1116 | +and likewise for $\beta_3$. The set $\mathcal{T}'$ contains the true event, so |
| |
1117 | +the smallest attainable value is $1/|\mathcal{T}'|$. This statistic is the |
| |
1118 | +``placebo in time'' diagnostic recommended for interrupted time series |
| |
1119 | +\citep{bernal2018controls}: if a level shift of comparable $t$-statistic is |
| |
1120 | +found at a quarter of all candidate dates, a conventional $p$-value from the |
| |
1121 | +HAC standard error of $\beta_2$ overstates the evidence, because the series |
| |
1122 | +has enough curvature and persistence that any date looks like a break. The |
| |
1123 | +value is not a $p$-value under a null hypothesis in the usual sense (the |
| |
1124 | +candidate dates are not exchangeable draws from a null distribution, and the |
| |
1125 | +regressions share the same data), and the manuscript calls it a placebo rank. |
| |
1126 | +Two implementation details depart from the plan's wording. The plan speaks of |
| |
1127 | +placebo breaks ``at every month of the baseline''; the implementation ranges |
| |
1128 | +over the baseline \emph{and} the post period, so that a change in 2023 or |
| |
1129 | +2024 that is larger than the one at the reference date lowers the rank of the |
| |
1130 | +reference date, which is the conservative direction for hypothesis H1 and is |
| |
1131 | +why the wider range was kept. And the minimum of 24 pre-months and 12 |
| |
1132 | +post-months means that candidate dates near the ends of the series, where the |
| |
1133 | +level-change coefficient is poorly identified, are excluded from the |
| |
1134 | +comparison. |
| |
1135 | + |
| |
1136 | +\subsection{Multiple breaks: dynamic programming with information criteria} |
| |
1137 | +\label{sec:m-bp} |
| |
1138 | + |
| |
1139 | +The Bai--Perron framework \citep{bai1998breaks, bai2003computation} estimates |
| |
1140 | +$m$ unknown break dates $\tau_1 < \dots < \tau_m$ as the global minimiser of |
| |
1141 | +the sum of squared residuals of a segment-wise regression. The implementation |
| |
1142 | +is the pure dynamic programme of \citet{bai2003computation} on the log series |
| |
1143 | +with a segment-wise linear trend: for a segment $[i, j)$ of at least $h$ |
| |
1144 | +observations, |
| |
1145 | +\begin{equation} |
| |
1146 | + S(i, j) \;=\; \min_{a, b}\ \sum_{t=i}^{j-1} \bigl( y_t - a - b\,t \bigr)^2, |
| |
1147 | +\end{equation} |
| |
1148 | +and the optimal partition of the first $j$ observations into $m+1$ segments |
| |
1149 | +is obtained recursively, |
| |
1150 | +\begin{equation} |
| |
1151 | + C_m(j) \;=\; \min_{\,m h \le i \le j - h}\ \bigl\{ C_{m-1}(i) + S(i, j) \bigr\}, \qquad C_0(j) = S(0, j), |
| |
1152 | + \label{eq:bp-dp} |
| |
1153 | +\end{equation} |
| |
1154 | +with the break positions recovered by back-tracking. The minimum segment |
| |
1155 | +length is $h = \max\{8, \lfloor 0.15\, n \rfloor\}$ months, a trimming of 15 per |
| |
1156 | +cent that is the usual choice in this literature; it forbids breaks in the |
| |
1157 | +first and last $h$ months, which for the default window means that no break |
| |
1158 | +can be placed after the autumn of 2024, and it also forbids two breaks closer |
| |
1159 | +than $h$ months. The number of breaks, up to three, is selected by two |
| |
1160 | +criteria computed on the global minima $\mathrm{SSR}_m = C_m(n)$, |
| |
1161 | +\begin{align} |
| |
1162 | + \mathrm{BIC}(m) &= n \log\bigl(\mathrm{SSR}_m / n\bigr) + p_m \log n, \\ |
| |
1163 | + \mathrm{LWZ}(m) &= n \log\bigl(\mathrm{SSR}_m / (n - p_m)\bigr) + p_m\, c_0\, (\log n)^{2.1}, \qquad c_0 = 0.299, |
| |
1164 | +\end{align} |
| |
1165 | +with $p_m = 2(m+1) + m$ parameters (an intercept and a slope per segment plus |
| |
1166 | +the $m$ break dates). The BIC is Schwarz's criterion \citep{schwarz1978bic}; |
| |
1167 | +the LWZ criterion of \citet{liu1997segmented}, in the form used by |
| |
1168 | +\citet{bai2003computation}, penalises additional breaks more heavily and is |
| |
1169 | +known to select fewer spurious breaks when the errors are serially correlated, |
| |
1170 | +which they are in these series (the monthly residual autocorrelation is what |
| |
1171 | +the HAC errors of the ITS correct for). Both selections are reported. The |
| |
1172 | +sequential sup-$F$ tests of \citeauthor{bai1998breaks} are not implemented; in |
| |
1173 | +their absence the break dates carry no $p$-value, and their stability is |
| |
1174 | +assessed instead through the agreement of the two criteria, the sensitivity to |
| |
1175 | +$t_{\mathrm{L}}$, and the agreement with the two methods below. The regression |
| |
1176 | +in each segment contains no seasonal terms, so a seasonal pattern can only be |
| |
1177 | +absorbed by the error; with a trimming of eight or more months per segment |
| |
1178 | +this cannot create a break by itself but it inflates $\mathrm{SSR}_m$ uniformly. |
| |
1179 | + |
| |
1180 | +\subsection{PELT} |
| |
1181 | + |
| |
1182 | +The pruned exact linear time algorithm \citep{killick2012pelt} finds the |
| |
1183 | +changepoint set $\mathcal{C}$ minimising a penalised cost, |
| |
1184 | +\begin{equation} |
| |
1185 | + \min_{\mathcal{C}}\ \sum_{\text{segments } [i,j) \text{ induced by } \mathcal{C}}\ \sum_{t=i}^{j-1} (y_t - \bar y_{i:j})^2 \;+\; \lambda\, |\mathcal{C}| , |
| |
1186 | + \qquad \lambda = 2\, \widehat{\mathrm{Var}}(y)\, \log n , |
| |
1187 | + \label{eq:pelt} |
| |
1188 | +\end{equation} |
| |
1189 | +with the $\ell_2$ cost (piecewise constant mean) and a minimum segment length |
| |
1190 | +of six months, as implemented in the \texttt{ruptures} library |
| |
1191 | +\citep{truong2020ruptures}. The penalty is of the Schwarz type studied by |
| |
1192 | +\citet{yao1988schwarz} for the piecewise-constant Gaussian model, scaled by |
| |
1193 | +the sample variance of the series; it is a conservative choice in a trending |
| |
1194 | +series, because a piecewise-constant mean must spend changepoints to follow a |
| |
1195 | +trend, so the reported PELT dates should be read as the points where the level |
| |
1196 | +moves most, not as the Bai--Perron break dates of a trend model. The two |
| |
1197 | +methods are complementary for that reason: a date flagged by both is a change |
| |
1198 | +in level that is not explained by a change in slope, and a date flagged only |
| |
1199 | +by the trend model is a kink. |
| |
1200 | + |
| |
1201 | +\subsection{Bayesian online changepoint detection} |
| |
1202 | +\label{sec:m-bocpd} |
| |
1203 | + |
| |
1204 | +The online detector of \citet{adams2007bocpd} maintains, after each |
| |
1205 | +observation, a posterior distribution over the \emph{run length} $r_t$, the |
| |
1206 | +number of months since the last changepoint. Let $H$ be the constant hazard |
| |
1207 | +(the prior probability that a change occurs at any month, here $1/60$, a prior |
| |
1208 | +mean run length of five years), and let the observations within a run be |
| |
1209 | +modelled as $y \sim \mathcal{N}(\mu, \sigma^2)$ with the conjugate |
| |
1210 | +Normal--Inverse-Gamma prior |
| |
1211 | +$\mu \mid \sigma^2 \sim \mathcal{N}(\mu_0, \sigma^2/\kappa_0)$, |
| |
1212 | +$\sigma^2 \sim \mathrm{IG}(\alpha_0, \beta_0)$. After $r$ observations |
| |
1213 | +$y_1, \dots, y_r$ of the current run the hyperparameters update as |
| |
1214 | +\begin{equation} |
| |
1215 | + \kappa_r = \kappa_0 + r, \quad |
| |
1216 | + \mu_r = \frac{\kappa_0 \mu_0 + \sum_{i \le r} y_i}{\kappa_r}, \quad |
| |
1217 | + \alpha_r = \alpha_0 + \frac{r}{2}, \quad |
| |
1218 | + \beta_r = \beta_0 + \tfrac12 \sum_{i \le r}(y_i - \bar y)^2 + \frac{\kappa_0\, r\, (\bar y - \mu_0)^2}{2\,\kappa_r}, |
| |
1219 | +\end{equation} |
| |
1220 | +which the code computes recursively, one observation at a time, as |
| |
1221 | +$\mu' = (\kappa\mu + y)/(\kappa+1)$, $\kappa' = \kappa + 1$, $\alpha' = \alpha + \tfrac12$, |
| |
1222 | +$\beta' = \beta + \kappa (y - \mu)^2 / (2(\kappa+1))$. The predictive |
| |
1223 | +distribution of the next observation given a run of length $r$ is the |
| |
1224 | +Student-$t$ |
| |
1225 | +\begin{equation} |
| |
1226 | + y_{t+1} \mid r_t = r \;\sim\; t_{2\alpha_r}\Bigl( \mu_r,\ \frac{\beta_r (\kappa_r + 1)}{\alpha_r \kappa_r} \Bigr), |
| |
1227 | + \label{eq:predictive} |
| |
1228 | +\end{equation} |
| |
1229 | +and the run-length recursion is |
| |
1230 | +\begin{align} |
| |
1231 | + P(r_{t+1} = r + 1,\ y_{1:t+1}) &= P(r_t = r,\ y_{1:t})\ p(y_{t+1} \mid r_t = r)\ (1 - H), \label{eq:rl-growth}\\ |
| |
1232 | + P(r_{t+1} = 0,\ y_{1:t+1}) &= \sum_{r} P(r_t = r,\ y_{1:t})\ p(y_{t+1} \mid r_t = r)\ H, \label{eq:rl-change} |
| |
1233 | +\end{align} |
| |
1234 | +normalised at each step to give the posterior $R_{t+1}(r) = P(r_{t+1} = r \mid y_{1:t+1})$. |
| |
1235 | +The hyperparameters are initialised from the first six months of the series |
| |
1236 | +($\mu_0$ their mean, $\beta_0$ their variance, $\kappa_0 = 1$, $\alpha_0 = 1$), |
| |
1237 | +an empirical-Bayes choice that uses the start of the series twice and whose |
| |
1238 | +effect vanishes after a few dozen observations. Two detection statistics are |
| |
1239 | +derived. Because a constant hazard places mass $H$ on $r = 0$ at every step by |
| |
1240 | +construction, the probability of a change ``at $t$'' is uninformative on its |
| |
1241 | +own; the code reports instead $P(r_t \le 3 \mid y_{1:t})$, the posterior |
| |
1242 | +probability that a change occurred within the last four months, and the |
| |
1243 | +maximum of this quantity over December 2022 to March 2023 is the statistic |
| |
1244 | +reported for the reference date. The second statistic is the sequence of |
| |
1245 | +maximum-a-posteriori run lengths $\hat r_t = \arg\max_r R_t(r)$; a |
| |
1246 | +\emph{reset} is declared at month $t$ when $\hat r_t \le 2$ after |
| |
1247 | +$\hat r_{t-1} > 12$, that is, when the most probable explanation of the data |
| |
1248 | +switches from a run of more than a year to a run that has just begun, and the |
| |
1249 | +list of reset months is the online detector's set of changepoints. The |
| |
1250 | +exact offline analogue \citep{fearnhead2006exact} would give the smoothed |
| |
1251 | +rather than the filtered posterior; the filtered version is used because the |
| |
1252 | +question it answers---when would a contemporaneous observer have detected the |
| |
1253 | +change---is the one relevant to the lead-time evaluation of |
| |
1254 | +Section~\ref{sec:m-emerging}. |
| |
1255 | + |
| |
1256 | +\subsection{Alternative event dates} |
| |
1257 | +\label{sec:m-altdates} |
| |
1258 | + |
| |
1259 | +The ITS is re-estimated for four alternative markers (June 2020 for GPT-3, |
| |
1260 | +May 2022 for DALL$\cdot$E~2, April 2023 for GPT-4, October 2024 for o1), each |
| |
1261 | +coded as the first full month after the public release, and an event is |
| |
1262 | +skipped when fewer than twelve months precede or follow it in the series. The |
| |
1263 | +exercise is not a horse race among causes. Its purpose is to show whether the |
| |
1264 | +estimated level and slope changes are specific to the November 2022 marker or |
| |
1265 | +whether the series exhibits similar changes at other dates associated with |
| |
1266 | +the same technological sequence, in which case the honest description is of a |
| |
1267 | +multi-year acceleration rather than a single break; the Bai--Perron dates are |
| |
1268 | +the data-driven complement to this pre-specified list. For the total series |
| |
1269 | +the sensitivity of the reference-date ITS and of the Bai--Perron dates to the |
| |
1270 | +truncation month $t_{\mathrm{L}}$ is also tabulated. |
| |
1271 | + |