index.Rmd

---
title: "Markdown - Breaking in or Breaking Through? How Local Specialisations Shape the Integration  of AI Technologies"
author: "Matheus E. Leusin"
date: "2024-09-12"
output: 
  html_document:
    toc: true
    toc_float: true # This makes it a floating sidebar in HTML
    theme: united
    df_print: paged
  word_document: default
  pdf_document: default
---

<br>

> 💡 **Welcome to the Methodological Appendix!**
>
> This file provides a fully transparent and reproducible account of the data analysis performed for the paper, *'Breaking in or Breaking Through? How Local Specialisations Shape the Integration  of AI Technologies'* (doi: https://doi.org/10.1080/10438599.2025.2558626). Its goal is to offer a clear roadmap of the methodological steps, from initial data sourcing to the final regression analysis.

<br>

This document details the complete analytical pipeline, structured into five key parts:

> **1. Calculating Technological Specializations** *(Section 1.1)* ⚙️<br>
><small> This section describes how raw patent data is processed to calculate **Revealed Technological Advantage (RTA)** for different countries and for AI as a whole. This is performed across three distinct time intervals to capture technological evolution.</small><br>
>**2. Constructing and Visualizing Technological Spaces** *(Sections 1.2 & 1.3)* 🌐<br>
><small> Here, we explain the creation of the **Global Technological Space (GTS)**, based on the co-occurrence of technological fields in patents, and the dynamic **AI-specific Technological Space (ATS)**. This section also covers the visualization of national trajectories within these spaces.</small><br>
>**3. Generating Supporting Figures** *(Section 2)* 📊<br>
><small> This part outlines the code used to create the descriptive figures presented in the paper, such as the growth of AI patents and the share of specialized fields.</small><br>
>**4. Robustness Checks & Permutation Analysis** *(Section 3)* 🎲<br>
><small> This section details the permutation-based robustness checks, which are used to ensure that our findings are statistically significant and not the result of random chance.</small><br>
>**5. Econometric Analysis** *(Section 4)* 📈<br>
><small> Finally, we present the regression models used to formally test our hypotheses, including the data setup, model specifications, and interpretation of the results.</small><br>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, error = FALSE)
```

# 1. Technological Spaces based on Technological field

## 1.1. Calculate Specializations for Different Time intervals

This section details the foundational step of our analysis: calculating the technological specializations of countries and of AI itself. We begin by loading extensive patent datasets and defining several custom functions to streamline the process. The end goal is to compute Revealed Technological Advantage (RTA) scores for three distinct time intervals, which will later serve as the basis for constructing our technological spaces.

**Data Loading and Initial Setup**

We start by loading the necessary R libraries and defining a set of custom functions that will be used repeatedly for data aggregation and weighting. The primary dataset is a large file containing patent applications and their associated inventor locations, which is loaded in manageable chunks to optimize memory usage.

```{r, include=FALSE}
#### Main Code
# Load required libraries
library(tidyverse)
library(magrittr)
library(tidygraph)
library(ggraph)
library(EconGeo)
library(data.table)
library(netrankr)
library(dplyr)
library(tidyr)
library(ggrepel)
library(scales)
library(patchwork)
library(RColorBrewer)
library(janitor)
library(ggforce)
library(stringr)
library(openxlsx)
library(kableExtra)

rm(list=ls())
#set the working directory to where you saved the R code:
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#1.First part: Technological Spaces based on Technological field -----
#Create important functions 
# This function groups data by application ID, and within each group, it calculates
# a field-specific weight equal to 1 divided by the number of records in that group.
group_by_applnID <- function(data){
  data %>%
    group_by(appln_id) %>%
    mutate(field_weight = 1 / n()) %>%
    ungroup()
}

# This function aggregates the weighted fields at the country-technology field level.
group_by_ctry_and_techn_field <- function(data){
  data %<>%
    group_by(ctry_code, techn_field_nr) %>%
    summarise(n_tech_reg = sum(field_weight)) %>%
    ungroup() %>%
    drop_na() 
}

# This function aggregates the weighted fields at the country-subclass level.
group_by_ctry_and_subclass <- function(data){
  data %<>%
    group_by(ctry_code, Subclass) %>%
    summarise(n_tech_reg = sum(field_weight)) %>%
    ungroup() %>%
    drop_na() 
}

# This function arranges multiple ggplots into one figure.
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Combine all plots into a single list
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  
  # If no layout is provided, define one based on the specified number of columns
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  # If only one plot, just print it
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up a new page for the layout
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Print each plot in the correct layout position
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

# This function creates a sparse matrix from the given inputs.
# i.input and j.input represent row and column indices, respectively.
create_sparse_matrix <- function(i.input, j.input){
  require(Matrix)
  mat <- spMatrix(
    nrow = i.input %>% n_distinct(),
    ncol = j.input %>% n_distinct(),
    i    = i.input %>% factor() %>% as.numeric(),
    j    = j.input %>% factor() %>% as.numeric(),
    x    = rep(1, i.input %>% length())
  )
  
  row.names(mat) <- i.input %>% factor() %>% levels()
  colnames(mat)  <- j.input %>% factor() %>% levels()
  return(mat)
}

rows_part2 <- 45182803 - 40000000
ipc_all_patents_part2_chunk1_df <- fread("large_files/All_patents_and_IPCs_Part2.csv", 
                                         header = FALSE, nrow = 20000000)
ipc_all_patents_part2_chunk2_df <- fread("large_files/All_patents_and_IPCs_Part2.csv", 
                                         header = FALSE, nrow = 20000000, skip = 20000000)
ipc_all_patents_part2_chunk3_df <- fread("large_files/All_patents_and_IPCs_Part2.csv", 
                                         header = FALSE, nrow = rows_part2, skip = 40000000)

ipc_all_patents_part2_df <- rbind(ipc_all_patents_part2_chunk1_df, 
                                  ipc_all_patents_part2_chunk2_df, ipc_all_patents_part2_chunk3_df)

rm(ipc_all_patents_part2_chunk1_df, ipc_all_patents_part2_chunk2_df, ipc_all_patents_part2_chunk3_df)

setnames(ipc_all_patents_part2_df, c("appln_id", "ctry_code", "techn_field_nr", "weight", "priority_year"))

ai_patents_df <- fread("other_files/IPCs_AI.csv", sep=";", dec=",", header=TRUE)
ai_patents_df$ctry_code <- "AI_pat"

# Load IPC technology names for labeling
ipc_names_df <- read.csv("other_files/ipc_technology.csv", sep = ";", header = TRUE) %>%
  select(field_nr, sector, field_name) %>% distinct(field_nr, .keep_all = TRUE) %>%
  mutate(techn_field_nr = field_nr) %>% arrange(techn_field_nr)
```

A preview of the primary patent data (`ipc_all_patents_part2_df`) is shown below. Each row represents a patent application (`appln_id`) linked to an inventor's country (`ctry_code`) and a technological field (`techn_field_nr`).

```{r}
kable(as.data.frame(ipc_all_patents_part2_df[1:6,]))
```
Note that the original `weight` column from PATSTAT is disregarded in our analysis. We recalculate a fractional weight internally to ensure that each patent application has a total weight of 1, distributed equally among its assigned technological fields.

Next, we load two supplementary datasets:

1. **AI Patent Data** (`other_files/IPCs_AI.csv`): A curated list of patent applications identified as being related to Artificial Intelligence.

2. **IPC Technology Names** (`other_files/ipc_technology.csv`): A reference file that maps technological field numbers to their descriptive names and sectors.

The *AI Patent Data* (`ai_patents_df`) looks like this:
```{r}
head(ai_patents_df)
```
The *IPC Technology Names* (`ipc_names_df`) looks like this:
```{r}
kable(as.data.frame(ipc_names_df[1:6,]))
```
**Calculating Specializations for Interval 1 (1974-1988)** 

The core of this section involves calculating the specialization scores for our first time interval, 1974-1988. This process is repeated identically for the subsequent two intervals.

The first step is to filter the main patent dataset for the specified period. We then apply our custom functions, `group_by_applnID()` and `group_by_ctry_and_techn_field()`, to fractionally count patent activities.

```{r, include=FALSE}
####1.1.2.1. Interval 1 (1974 - 1988)
# Define Time Interval
start_year <- 1973
end_year   <- 1989

###Filter Data for Interval 1 
ipc_all_patents_first_period_df <- ipc_all_patents_part2_df[priority_year > start_year & priority_year < end_year]
ipc_all_patents_first_period_df <- ipc_all_patents_first_period_df[, .(appln_id, ctry_code, techn_field_nr)]

# Calculate General RCA for Interval 1
# Compute weighted values at country-technology field level
region_tech_fields_1_df <- group_by_applnID(ipc_all_patents_first_period_df)
```

The `group_by_applnID()` function assigns an equal weight to each technological field within a single patent. For instance, if a patent is classified under four fields, each field receives a weight of 0.25. The result is a weighted dataset:

```{r}
kable(as.data.frame(region_tech_fields_1_df[1:6,]))
```

Next, `group_by_ctry_and_techn_field()` aggregates these weights, summing them up for each country-technology pair. This yields the total fractional count of patents for each country in each technological field.

```{r, include=FALSE}
region_tech_fields_1_df <- group_by_ctry_and_techn_field(region_tech_fields_1_df)
```

```{r}
kable(as.data.frame(region_tech_fields_1_df[1:6,]))
```

This aggregated data, saved as `reg_tech_FirstPeriod.csv`, is transformed into a country-technology matrix. The matrix rows represent countries, columns represent technological fields, and the values are the fractional patent counts. It looks like this:

```{r, include=FALSE}
mat_reg_tech1 <- region_tech_fields_1_df %>% arrange(techn_field_nr, ctry_code) %>%
  pivot_wider(names_from = techn_field_nr, values_from = n_tech_reg, values_fill = 0)
mat_reg_tech1 %<>% remove_rownames %>% column_to_rownames(var="ctry_code") %>%
  as.matrix() %>% round()
```

```{r}
kable(as.matrix(mat_reg_tech1[1:20, 1:12]), caption = "Sample of the Country-technology matrix")
```

Finally, we use this matrix to calculate the **Revealed Technological Advantage (RTA)** for each country in each field. RTA is a non-binary index that measures whether a country has a greater share of patents in a specific technology compared to the global average. An RTA value greater than or equal to 1 indicates a specialization.

```{r, include=FALSE}
## Calculate RCA (relative comparative advantage) for general technologies
reg_RCA1_df <- mat_reg_tech1 %>% location_quotient(binary = FALSE) %>% as.data.frame() %>% 
  rownames_to_column("ctry_code") %>% as_tibble() %>% gather("techn_field_nr", "RCA", -ctry_code)
```

```{r}
kable(as.data.frame(reg_RCA1_df[1:6,]))
```

This entire process is then repeated, this time using only the AI-related patents from the first interval to calculate **AI-specific RTAs** for each country. For the AI patents, the RTAs look like this:

```{r, include=FALSE}
# Calculate AI-Specific RCA for Interval 1
ai_patents_period_1_df <- ai_patents_df[ai_patents_df$priority_year > start_year & ai_patents_df$priority_year < end_year,]

# Join technological fields info to AI patents
ai_patents_period_1_df <- distinct(ai_patents_period_1_df, appln_id, .keep_all = TRUE)[, c(1,3)]
ai_patents_period_1_df <- left_join(ai_patents_period_1_df, ipc_all_patents_first_period_df, by = "appln_id")

# Compute weighted fields for AI patents
region_tech_fields_1_ai_df <- group_by_applnID(ai_patents_period_1_df)
rm(ai_patents_period_1_df)
region_tech_fields_1_ai_df <- group_by_ctry_and_techn_field(region_tech_fields_1_ai_df)

mat_reg_tech1_AI <- region_tech_fields_1_ai_df %>%
  arrange(techn_field_nr, ctry_code) %>%
  pivot_wider(names_from = techn_field_nr, values_from = n_tech_reg, values_fill = list(n_tech_reg = 0))

mat_reg_tech1_AI %<>% remove_rownames %>% column_to_rownames(var="ctry_code") %>% as.matrix() %>%  round()

reg_RCA1_AI_df <- mat_reg_tech1_AI %>% location_quotient(binary = FALSE) %>% 
  as.data.frame() %>% rownames_to_column("ctry_code") %>%   
  as_tibble() %>%   gather(key = "techn_field_nr", value = "RCA", -ctry_code) %>%  arrange(ctry_code, techn_field_nr)
```

```{r}
kable(as.data.frame(reg_RCA1_AI_df[1:12,]))
```

The general and AI-specific RTA dataframes are then merged. The resulting file for the first interval shows, for each country and technological field, both its general specialization (`RCA_Gen`) and its AI-specific specialization (`RCA_AI`), as highlighted below for the whole dataset, and for Japan as an example.
```{r, include=FALSE}
# Merge general and AI-specific RCA data for Interval 1
rca_data_period_1_df <- merge(reg_RCA1_df, reg_RCA1_AI_df, all = TRUE, by = c("ctry_code", "techn_field_nr"))
rca_data_period_1_df$Period <- "1974-1988"
names(rca_data_period_1_df) <- c("ctry_code", "techn_field_nr", "RCA_Gen", "RCA_AI", "Period")
```

```{r}
#Resulting file:
kable(as.data.frame(rca_data_period_1_df[1:6,]))
#Example Japan:
kable(as.data.frame(rca_data_period_1_df[rca_data_period_1_df$ctry_code == "JP",][1:6,]))
```

A key methodological step follows: we treat the entire corpus of AI patents as if it belonged to a single, hypothetical 'country' named `AI_pat`. This novel approach allows us to calculate the RTA for AI itself across all technological fields, providing a benchmark against which national specializations can be compared. The resulting data is saved for later use, and it looks like this:

```{r, include=FALSE}
# Regional Tech for AI - Interval 1
setDT(ai_patents_df)
setDT(ipc_all_patents_first_period_df)
ipc_all_patents_first_period_df[ai_patents_df, on = "appln_id", ctry_code := i.ctry_code]
region_tech_ai_1_df <- group_by_applnID(ipc_all_patents_first_period_df)
region_tech_ai_1_df <- group_by_ctry_and_techn_field(region_tech_ai_1_df)
```

```{r}
kable(as.data.frame(region_tech_ai_1_df[region_tech_ai_1_df$ctry_code == "AI_pat",][1:6,]))
```

**Consolidating Data Across All Intervals**

The calculation process detailed above is repeated for the remaining two intervals: **1989-2003** and **2004-2018**. After processing all periods, the three interval-specific RTA files are combined into a single, comprehensive dataset named `IPC_RCAs` that is saved for later usage  (`Files_created_with_the_code/data/files_code_Fields_analysis/IPC_RCAs.csv`). This file contains the complete history of general and AI-specific specializations for all countries across the three time periods. 
Using Japan again as an example, the file looks like this for this country across each interval:

```{r, include=FALSE}
# Clean up
IPC_RCAs <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/IPC_RCAs.csv", 
                     sep = ";", header = TRUE, dec=",")
```

```{r}
kable(as.data.frame(IPC_RCAs[IPC_RCAs$ctry_code == "JP" & IPC_RCAs$Period == "1974-1988",][1:6,]))
kable(as.data.frame(IPC_RCAs[IPC_RCAs$ctry_code == "JP" & IPC_RCAs$Period == "1989-2003",][1:6,]))
kable(as.data.frame(IPC_RCAs[IPC_RCAs$ctry_code == "JP" & IPC_RCAs$Period == "2004-2018",][1:6,]))
```

To facilitate further analysis and visualization, we create consolidated summary files for each interval. These files combine the RTA scores of the four focus countries (US, CN, KR, JP) and the `AI_pat` entity into a single, wide-format table. They are named, e.g., as `Files_created_with_the_code/data/files_code_Fields_analysis/Metrics_First_period.csv` (for the First interval, and the names change accordingly for the Second and Third intervals). The data for this example interval looks like this:

```{r, include=FALSE}
###1.1.2.4.  Create a Specializations Summary
# First Interval Countries 
# Load first interval data for countries
reg_tech1_countries <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/reg_tech_FirstPeriod.csv", 
                                sep = ";", header = TRUE, dec=",")

# Create a wide matrix of technology fields for the first interval
mat_reg_tech1_countries <- reg_tech1_countries %>%  arrange(techn_field_nr, ctry_code) %>%
  pivot_wider(names_from = techn_field_nr, values_from = n_tech_reg, values_fill = list(n_tech_reg = 0))

mat_reg_tech1_countries %<>%   remove_rownames %>% 
  column_to_rownames(var="ctry_code") %>%  as.matrix() %>%    round()

# Compute RCA for the first interval (general)
reg_RCA1_countries <- mat_reg_tech1_countries %>%   location_quotient(binary = FALSE) %>% 
  as.data.frame() %>% rownames_to_column("ctry_code") %>% 
  as_tibble() %>% gather(key = "techn_field_nr", value = "RCA", -ctry_code) %>% arrange(ctry_code, techn_field_nr)

# First Interval AI 
# Load first interval data for AI-specific patents
reg_tech1_AI <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/reg_techAI_FirstPeriod.csv", 
                         sep = ";", header = TRUE, dec=",")

# Create a wide matrix of technology fields for AI patents in the first interval
mat_reg_tech1_AI <- reg_tech1_AI %>%  arrange(techn_field_nr, ctry_code) %>%
  pivot_wider(names_from = techn_field_nr, values_from = n_tech_reg, values_fill = list(n_tech_reg = 0))

mat_reg_tech1_AI %<>%   remove_rownames %>%   column_to_rownames(var="ctry_code") %>%
  as.matrix() %>%  round()

# Compute RCA for the first interval (AI-specific)
reg_RCA1_AI <- mat_reg_tech1_AI %>%   location_quotient(binary = FALSE) %>% 
  as.data.frame() %>%   rownames_to_column("ctry_code") %>% 
  as_tibble() %>%   gather(key = "techn_field_nr", value = "RCA", -ctry_code) %>%
  arrange(ctry_code, techn_field_nr)

# Load IPC names and metadata
IPC_names <- read.csv("other_files/ipc_technology.csv", sep = ";", header = TRUE)%>%
  select(field_nr, sector, field_name) %>%  distinct(field_nr, .keep_all = TRUE) %>%
  mutate(techn_field_nr = field_nr) %>%  arrange(techn_field_nr)
IPC_names <- IPC_names[, -1]

# Extract top countries and AI data for first interval
US_first_period <- reg_RCA1_countries[,2:3][reg_RCA1_countries$ctry_code == "US",]
CN_first_period <- reg_RCA1_countries[,2:3][reg_RCA1_countries$ctry_code == "CN",]
KR_first_period <- reg_RCA1_countries[,2:3][reg_RCA1_countries$ctry_code == "KR",]
JP_first_period <- reg_RCA1_countries[,2:3][reg_RCA1_countries$ctry_code == "JP",]
AI_first_period <- reg_RCA1_AI[,2:3][reg_RCA1_AI$ctry_code == "AI_pat",]

# Merge IPC names with the countries and AI RCAs for the first interval
First_period <- merge(merge(merge(merge(merge(IPC_names, US_first_period), CN_first_period, by = "techn_field_nr"), 
  KR_first_period, by = "techn_field_nr"), JP_first_period, by = "techn_field_nr"), 
  AI_first_period, by = "techn_field_nr")

names(First_period) <- c("techn_field_nr", "sector", "field_name", "RCA_US", "RCA_CN","RCA_KR","RCA_JP", "RCA_AI")
```

```{r}
kable(as.data.frame(First_period[1:6,]))
```

Finally, these three interval-specific summary files are merged into a master file named `All_periods`, shown below. This file includes additional labels for analytical purposes, though these are not central to the paper's main findings.

```{r, include=FALSE}
IPC_names <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Specializations_All_periods_IPC.csv", 
                      sep = ";", header = TRUE)[,c(-1)]
```

```{r}
head(IPC_names)
```

In the last step of this sub-section, we use the `IPC_RCAs.csv` file to generate a summary table (`IPC_RCAs_Top4`). Here, the non-binary RTA values are binarized, where any RTA ≥ 1 is considered a specialization (value of 1) and any RTA < 1 is not (value of 0). We then sum these binary indicators to count the number of general specializations, AI-specific specializations, and coinciding specializations (where a country is specialized in both the general field and its AI-specific application) for each country and interval, resulting in the following dataset:

```{r, include=FALSE}
IPC_RCAs_Top4 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/RCA_4countries_detailed.csv", 
                          sep = ";", header = TRUE, dec=",")
```

```{r}
kable(as.data.frame(IPC_RCAs_Top4[1:6,]))
```

## 1.2. Building the Global Technological Space (GTS)

The next step is to construct the backbone of our analysis: the **Global Technological Space (GTS)**. This space is a network where nodes represent technological fields, and the links between them signify their relatedness. We measure this relatedness based on the principle that technologies that frequently appear together within the same patent are likely to be related.

### 1.2.1. From Patents to a Co-occurrence Matrix

To quantify this relationship, we must first count how often every possible pair of technologies co-occurs across the entire patent dataset. We start by loading the complete patent database (which, due to its size, is again handled in chunks) and applying the `create_sparse_matrix` function. This function generates a very large matrix where rows are unique patents and columns are the 35 technological fields.

```{r, include=FALSE}
c <- 58841893 - 40000000
IPC_all_patents_Part1 <- fread("large_files/All_patents_and_IPCs_Part1.csv", header = FALSE, nrow = 20000000)
IPC_all_patents_Part2 <- fread("large_files/All_patents_and_IPCs_Part1.csv", header = FALSE, nrow = 20000000, skip = 20000000)
IPC_all_patents_Part3 <- fread("large_files/All_patents_and_IPCs_Part1.csv", header = FALSE, nrow = c, skip = 40000000)

names(IPC_all_patents_Part1) <- c("appln_id", "ctry_code", "techn_field_nr", "weight", "priority_year")
names(IPC_all_patents_Part2) <- c("appln_id", "ctry_code", "techn_field_nr", "weight", "priority_year")
names(IPC_all_patents_Part3) <- c("appln_id", "ctry_code", "techn_field_nr", "weight", "priority_year")

# Create sparse matrices and compute cross-products for the first big file
mat_tech_AI1 <- create_sparse_matrix(i = IPC_all_patents_Part1 %>% pull(appln_id),
                                     j = IPC_all_patents_Part1 %>% pull(techn_field_nr)) #
```
The resulting sparse matrix, `mat_tech_AI1`, indicates the presence of a technology in a given patent, and it looks like this:
```{r}
kable(as.matrix(mat_tech_AI1[1:20, 1:12]), caption = "Sample of the Sparse AI matrix")
```

By calculating the cross-product of this matrix (`t(M) %*% M`), we transform it into a 35x35 square **co-occurrence matrix**. Each cell (`i, j`) in this new matrix contains a count of how many patents simultaneously list technology `i` and technology `j`. This square matrix looks like this:
```{r, include=FALSE}
mat_tech_AI1 %<>%   crossprod() %>%  as.matrix() 
```

```{r}
kable(as.matrix(mat_tech_AI1[1:35, 1:35]), caption = "Sample of the co-occurrence matrix")
```

After processing all data chunks, the individual co-occurrence matrices are summed to create a final, comprehensive matrix, which is then saved as `Matrix_IPC.csv`. This file looks like this:

```{r, include=FALSE}
matrix2 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Matrix_IPC.csv", 
                    sep = ";", header = FALSE)

matrix2 <- matrix2 %>%  row_to_names(row_number = 1)
matrix <- matrix2[,-1]
rownames(matrix) <- matrix2[,1]
matrix <- as.matrix(matrix)
mat_tech_AI_Final <- matrix
```

```{r}
kable(as.matrix(mat_tech_AI_Final[1:35, 1:35]))
```

### 1.2.2. Calculating Relatedness and Defining the Network

Raw co-occurrence counts can be misleading, as highly prevalent technologies will naturally co-occur more often with others, inflating their apparent relatedness. To correct for this, we normalize the matrix using the `relatedness()` function from the `EconGeo` package, which employs a cosine similarity index. The result is a relatedness matrix, where each value represents the strength of the relationship between two technologies. It looks like this:

```{r, include=FALSE}
# Calculate relatedness using cosine similarity
mat_tech_rel_AI <- mat_tech_AI_Final %>% relatedness(method = "cosine")
```

```{r}
kable(as.matrix(mat_tech_rel_AI[1:35, 1:35]))
```

With the relatedness matrix complete, we can now treat it as an adjacency matrix to build a network graph (`g_tech_AI`). The nodes' centrality (**Eigenvector centrality**) is calculated to determine their importance in the network. For visual clarity in later plots, links with below-average weight (relatedness) are filtered out. Finally, a Fruchterman-Reingold layout algorithm is applied to determine the spatial coordinates (`coords_tech_AI`) of each node for visualization, which results in the following coordinates:

```{r, include=FALSE}
# Load IPC names and categories
IPC_names <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Specializations_All_periods_IPC.csv", 
                      sep = ";", header = TRUE)%>%
  select(techn_field_nr, sector, field_name, Category) %>%  distinct(techn_field_nr, .keep_all = TRUE) %>%
  mutate(techn_field_nr = techn_field_nr) %>%  arrange(techn_field_nr)

# Build a graph from the relatedness matrix
g_tech_AI <- mat_tech_rel_AI %>%   as_tbl_graph(directed = FALSE) %N>%
  left_join(IPC_names %>% mutate(techn_field_nr = as.character(techn_field_nr)), 
            by = c("name" = "techn_field_nr")) %>%  mutate(dgr = centrality_eigen(weights = weight)) %E>%
  filter(weight >= mean(weight))

# Layout for visualization (Fruchterman-Reingold)
coords_tech_AI <- g_tech_AI %>%   igraph::layout.fruchterman.reingold() %>% as_tibble()
colnames(coords_tech_AI) <- c("x", "y")
# Alternatively, load predefined coordinates
coords_tech_AI <- read.csv("other_files/coords_tech_AI_layout1.csv", sep = ";", header = TRUE, dec=",")
```

```{r}
kable(as.data.frame(coords_tech_AI[1:10,]))
```

### 1.2.3. Preparing Data for Visualization

In the final step of this section, we prepare the specialization data (calculated in Section 1.1) for plotting onto the GTS. We load the summary file (`RCA_4countries_detailed.csv`) and create a new categorical variable (`Var1`) that classifies each country-technology pair into one of four states: no specialization (0), general specialization (1), AI-specific (break-through) specialization (2), or coinciding (break-in) specialization (3). This will allow us to map the countries' technological trajectories directly onto the GTS structure in the next section. The dataset looks like this:

```{r, include=FALSE}
# Load top 4 countries RCA details
IPC_RCAs_Top4 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/RCA_4countries_detailed.csv", 
                          sep = ";", header = TRUE, dec=",")
IPC_RCAs_Top4$Total_RCA <- as.factor(IPC_RCAs_Top4$Total_RCA)
IPC_RCAs_Top4$Period_sim <- as.numeric(factor(IPC_RCAs_Top4$Period, levels=unique(IPC_RCAs_Top4$Period)))
IPC_RCAs_Top4$techn_field_nr <- as.character(IPC_RCAs_Top4$techn_field_nr)

# Load AI RCA data and merge with IPC_RCAs_Top4
AI_RCA <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Specializations_All_periods_IPC.csv", 
                   sep = ";", header = TRUE, dec=",")
AI_RCA$Period_sim <- as.numeric(factor(AI_RCA$Period, levels=unique(AI_RCA$Period)))
AI_RCA <- AI_RCA[, c(2,9,13)]
AI_RCA$techn_field_nr <- as.character(AI_RCA$techn_field_nr)
names(AI_RCA) <- c("techn_field_nr", "RCA_AI_Period", "Period_sim")

IPC_RCAs_Top4 <- left_join(IPC_RCAs_Top4, AI_RCA, by = c("techn_field_nr", "Period_sim"))

# Adjust Total_RCA_2 to differentiate between general and AI specialization
IPC_RCAs_Top4$Total_RCA_2 <- IPC_RCAs_Top4$Round_general + 2*IPC_RCAs_Top4$Round_AI

# Summarize data by category of specialization
Newtable <- as.data.frame(table(IPC_RCAs_Top4$Total_RCA_2, IPC_RCAs_Top4$ctry_code, IPC_RCAs_Top4$Period))
Newtable$Var1 <- gsub("0", "No specialization", str_trim(Newtable$Var1))
Newtable$Var1 <- gsub("1", "General specialization", str_trim(Newtable$Var1))
Newtable$Var1 <- gsub("2", "AI-specific specialization", str_trim(Newtable$Var1))
Newtable$Var1 <- gsub("3", "Coinciding specialization", str_trim(Newtable$Var1))
```

```{r}
kable(as.data.frame(Newtable[1:10,]))
```

## 1.3. Plotting technological spaces

Now that the underlying data and network structures are in place, this section focuses on their visualization. We will generate the key plots presented in the paper, illustrating both the static, global structure of technology and the dynamic, evolving space of AI.

### 1.3.1. Global technological space (GTS)

We begin by plotting the fundamental structure of the Global Technological Space. This initial visualization is **geography-agnostic**, meaning it shows the inherent relatedness between technological fields without any country-specific data. The node size corresponds to its centrality (degree), and nodes are clustered and colored by their broader technological sector. This plot serves as the canvas upon which we will later map national trajectories.


```{r, fig.width=14, fig.height=10}
  g_tech_AI %>%  ggraph(layout =  coords_tech_AI) + 
  geom_edge_link(aes(width = weight), alpha = 0.4, colour = "grey") + 
  geom_node_point(aes(fill = sector, size = 1000^dgr, shape= sector))+ # 
  scale_shape_manual(values=c(21, 22, 23, 24, 25)) + scale_size("Degree", range = c(2, 12)) + 
  geom_node_text(aes(label = paste0(field_name, "\n(", name, ")")), size = 4, repel = TRUE) +  #field_name or name
  theme_graph(base_family = "sans")+  ggtitle("Global technological space: IPC Technological fields") + 
  theme(legend.title = element_text(size = 14), legend.text = element_text(size = 10)) + 
  guides(colour = guide_legend(override.aes = list(size=10)))+
  geom_mark_hull(aes(x = x, y=y, colour = sector, fill= sector,
                     linetype = sector), alpha = 0.15, expand = unit(2.5, "mm"), size = 1) 

```

Next, we overlay the country-specific specialization data onto the static GTS canvas. This allows us to visualize the technological trajectory of each country over the three time intervals. The shape of each node indicates the type of specialization (general, break-through, or break-in), while hulls are drawn to highlight the clusters of specialization for each period. This composite visualization reveals how each nation's technological focus has evolved within the global structure. Additionally, an horizontal bar plot is also generated to summarize the main indicators based on the country-specific specializations. Picking China as an example, this overlaid visualization and its linked bar-plot look like this:

```{r, fig.width=14, fig.height=10}
#GTS with specialisations per country
country_select <- c("CN", "US", "JP", "KR")
### 1.2.3.3. Third Country
i=1
IPC_RCAs_wide_simplified <- IPC_RCAs_Top4 %>% pivot_wider(id_cols = c(ctry_code, techn_field_nr, Label), 
    names_from = Period_sim,
    values_from = c(RCA_AI_Period, Total_RCA_2, RCA_Gen, RCA_AI, Round_general, Round_AI, Total_RCA), 
    names_glue = "{.value}_Period_{Period_sim}" )

  g_tech_AI %N>% left_join(IPC_RCAs_wide_simplified %>%
                             filter(ctry_code == country_select[i]) %>%
                             select(-ctry_code), by = c("name" = "techn_field_nr")) %>%
  mutate(Shape_Group_P1_Factor = factor(
    ifelse(is.na(Total_RCA_2_Period_1), "NA_Value", as.character(Total_RCA_2_Period_1)),
    levels = c("0", "1", "2", "3", "NA_Value"))) %>% ggraph(layout = coords_tech_AI) +
  geom_edge_link(aes(width = weight), alpha = 0.2, colour = "#CCCCCC", show.legend = FALSE) + 
  geom_node_point(aes(shape = Shape_Group_P1_Factor, 
                      size = 5, stroke = ifelse(Total_RCA_2_Period_1 == 3, 2.5, 1.3),
                      alpha = 1), color = "#FF3300", show.legend = c(shape=TRUE, size=FALSE, stroke=FALSE, alpha=FALSE, color=FALSE)) + 
  geom_node_point(aes(shape = factor(Total_RCA_2_Period_2),
                      size = 5.5, stroke = ifelse(Total_RCA_2_Period_2 == 3, 2.5, 1.3),
                      alpha = 1), color = "#3399FF", show.legend = FALSE) +
  geom_node_point(aes(shape = factor(Total_RCA_2_Period_3), 
                      size = 6.5,stroke = ifelse(Total_RCA_2_Period_3 == 3, 2.5, 1.3),
                      alpha = 1), color = "#009900", show.legend = FALSE) +
  scale_shape_manual(name = "Type of specialisation",
                     values = c("0" = 4, "1" = 1, "2" = 5, "3" = 2, "NA_Value" = 16), breaks = c("0", "1", "2", "3"),                                
                     labels = c("0" = "No specialisation", "1" = "General specialisation", 
                                "2" = "Break-through specialisation", "3" = "Break-in specialisation"), 
                     na.translate = FALSE, drop = FALSE) + scale_size("Degree", range = c(7, 18))+ 
  scale_alpha(guide = "none") + 
  #geom_node_label(aes(label = name), size = 2, repel = F) + 
  geom_mark_hull(aes(filter = Total_RCA_2_Period_1 > .99, x = x, y = y, fill = "Period 1", group = "Period 1"), 
                 concavity = .1, alpha = .11, linetype = "dotted",expand = unit(2, "mm"), size = .5, color = "#FF3300") + 
  geom_mark_hull(aes(filter = Total_RCA_2_Period_2 > .99, x = x, y = y, fill = "Period 2", group = "Period 2"),
                 concavity = .1, alpha = .11, linetype = "longdash",expand = unit(2, "mm"), size = .5, color = "#3399FF") +
  geom_mark_hull(aes(filter = Total_RCA_2_Period_3 > .99, x = x, y = y, fill = "Period 3", group = "Period 3"),
                 concavity = .1, alpha = .02, expand = unit(2, "mm"), size = 1, color = "#009900") +
  scale_fill_manual(name = "Interval colour (same for \nboth nodes and cluster)", # New legend for fill
                    values = c("Period 1" = "#FF3300", "Period 2" = "#3399FF", "Period 3" = "#009900"),
                    labels = c("Interval 1 (1974-1988)", "Interval 2 (1989-2003)", "Interval 3 (2004-2018)")) +
  theme_graph(base_family = "sans") +  theme(legend.position = "bottom", #right
                                             legend.box = "vertical", legend.title = element_text(size = 12, face = "bold"), 
                                             legend.text = element_text(size = 10), legend.key.size = unit(0.7, "cm") ) +
  ggtitle("d) Global technological space: China (1974-2018)") +
  geom_node_text(aes(label = name), size = 5, repel = TRUE) +  #field_name or name
  guides(shape = guide_legend(title.position = "top", 
                              override.aes = list(size = 5, stroke = 1.5, color = "black") ),
         colour = guide_legend(title.position = "top", 
                               override.aes = list(linetype = c("solid", "longdash", "dotted"), 
                                                   alpha = 1, size = 1, shape = NA) ))
  
  bar_plot_China <- bar_plot_China <- IPC_RCAs_Top4[IPC_RCAs_Top4$ctry_code == country_select[i],] %>%                                   
  arrange(Label, Period) %>%  group_by(Label) %>%                          
  mutate( general = Total_RCA_2 == 1,    
          break_in              = Total_RCA_2 == 2,            
          break_through         = Total_RCA_2 == 3,            
          sustained_general    = general  & lag(general, 1, default = FALSE),
          sustained_break_in    = break_in  & lag(break_in, 1, default = FALSE),
          sustained_break_through    = break_through & lag(break_through, 1, default = FALSE)) %>% 
  ungroup()

bar_plot_China <- bar_plot_China %>% 
  group_by(Period) %>% summarise(`General case`                 = sum(general,           na.rm = TRUE),
                                 `Break-through case`                 = sum(break_in,           na.rm = TRUE),
                                 `Break-in case`            = sum(break_through,      na.rm = TRUE),
                                 `Sustained General case`       = sum(sustained_general, na.rm = TRUE),
                                 `Sustained break-through case`       = sum(sustained_break_in, na.rm = TRUE),
                                 `Sustained break-in case`  = sum(sustained_break_through, na.rm = TRUE),
                                 .groups = "drop") %>% arrange(Period)

plot_long_China <- bar_plot_China |>  rename(Period = Period) |>
  pivot_longer(cols= -Period,names_to= "Indicator",values_to = "Count")

#order labels
plot_long_China$Indicator <- factor(plot_long_China$Indicator, levels = rev(c("General case", "Break-through case",  "Break-in case", 
                                                              "Sustained General case", "Sustained break-through case", "Sustained break-in case")))
plot_long_China$Period <- factor(plot_long_China$Period, levels = c("2004-2018", "1989-2003", "1974-1988"))

legend_order <- c(
  "General case", "Break-through case", "Break-in case",
  "Sustained General case", "Sustained break-through case", "Sustained break-in case"
)

  ggplot(plot_long_China, aes(x = factor(Period),y = Count, fill = Indicator)) +
  geom_col(position = position_dodge(width = .8), width = .7) +
  scale_fill_manual(values = c("General case"  = "#FF3300",
    "Sustained General case"  = "#993333",
    "Break-in case"                 = "#009900", #3399FF
    "Sustained break-in case"       = "#006633", #3333CC
    "Break-through case"            = "#3399FF",  #009900
    "Sustained break-through case"  = "#3333CC"),
    breaks = legend_order) + #006633
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(x = "Interval",y = "Number of cases", fill = NULL, title = NULL)+
  ggtitle("Summary of specialisations China") +
  theme_classic(base_size = 11) + theme(legend.position = "bottom")+ coord_flip()
```

The plotting code is structured to iterate through each of the four focus countries by changing the `i` variable. The resulting figures, each depicting a single country's trajectory over three periods alongside a summary bar chart, are then saved.

### 1.3.2. AI-specific technological space (ATS)

Unlike the static GTS, the **AI-specific Technological Space (ATS)** is dynamic. Its structure is recalculated for each time interval, reflecting the rapid evolution of AI technology. Here, the relatedness between fields is based only on their co-occurrence within AI patents for that specific period. This approach allows us to observe which technological fields form the core of AI innovation at different points in time.

```{r, include=FALSE}
rm(list = ls()[!sapply(ls(), function(x) is.function(get(x)))])
gc()
# ATS First interval
patents_AI_specific_1st <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/AI_ALL_patents.csv", sep = ";", header = TRUE, dec=",")

a = 1973
b = 1989

patents_AI_specific_1st <- patents_AI_specific_1st[patents_AI_specific_1st$priority_year < b,]
patents_AI_specific_1st <- patents_AI_specific_1st[patents_AI_specific_1st$priority_year > a,]

length(unique(patents_AI_specific_1st$appln_id)) #436
patents_AI_specific_1st <- patents_AI_specific_1st[is.na(patents_AI_specific_1st$appln_id)==F,]

mat_tech_AI <- create_sparse_matrix(i = patents_AI_specific_1st %>% pull(appln_id),
                                    j = patents_AI_specific_1st %>% pull(techn_field_nr))

mat_tech_AI %<>%   crossprod() %>%   as.matrix() 

mat_tech_rel_AI <- mat_tech_AI %>%   relatedness(method = "cosine")

IPC_names <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Specializations_All_periods_IPC.csv", sep = ";", header = TRUE)%>%
  select(techn_field_nr, sector, field_name, Category) %>%  distinct(techn_field_nr, .keep_all = TRUE) %>%
  mutate(techn_field_nr = techn_field_nr) %>%  arrange(techn_field_nr)

g_tech_AI <- mat_tech_rel_AI %>% as_tbl_graph(directed = FALSE) %N>%
  left_join(IPC_names %>% mutate(techn_field_nr = techn_field_nr %>% as.character()), by = c("name" = "techn_field_nr")) %>%
  mutate(dgr = centrality_eigen(weights = weight)) %E>%  filter(weight >= mean(weight))

#Create the Coordinates
coords_tech_AI <- g_tech_AI %>% igraph::layout.fruchterman.reingold() %>% as_tibble()
colnames(coords_tech_AI) <- c("x", "y")
```
Starting with the first interval (1974-1988), the top 10 most central technological fields in the AI space are:
```{r}
g_tech_AI %N>%   arrange(desc(dgr)) %>%  as_tibble() %>%  slice(1:10)
```

We use the previously calculated AI specialization data (`AI_RCA`) to highlight the core technologies in each period. A binary flag indicates whether AI has an RTA ≥ 1 in a given field (and Period_sim refers to each interval, going from 1 to 3), like this:

```{r, include=FALSE}
AI_RCA <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Specializations_All_periods_IPC.csv", sep = ";", header = TRUE, dec=",")
AI_RCA$Period_sim <- as.numeric(factor(AI_RCA$Period,levels=unique(AI_RCA$Period)))
AI_RCA <- AI_RCA[,c(2,9,13)]
AI_RCA$techn_field_nr <- as.character(AI_RCA$techn_field_nr)
names(AI_RCA) <- c("techn_field_nr", "RCA_AI_Period", "Period_sim")
AI_RCA$Binary <- ifelse(AI_RCA$RCA_AI_Period < 1, 0,1)
```

```{r}
kable(as.data.frame(AI_RCA[1:6,]))
```

The following code generates the ATS for the first interval (1974-1988). The nodes with labels are those where AI is specialized (RTA ≥ 1).

```{r, fig.width=14, fig.height=10}
AI_RCA1 <- AI_RCA[AI_RCA$Period_sim == 1,]
p=1
  g_tech_AI %N>%
  left_join(AI_RCA1 %>% filter(Period_sim == p), by = c("name" = "techn_field_nr")) %>%
  ggraph(layout = coords_tech_AI) + 
  geom_edge_link(aes(width = weight), alpha = 0.2, colour = "#CCCCCC") +
  geom_node_point(aes(fill = sector, size = 1000^dgr, shape= sector)) +
  scale_shape_manual(values=c(21, 22, 23, 24, 25)) + labs(color   = "RCA")+ scale_size("Degree", range = c(2, 12)) +
  geom_node_text(aes(filter=Binary > .99, label = field_name), size = 6, repel = TRUE) +
  theme_graph(base_family = "sans") + guides(colour = guide_legend(override.aes = list(size=5)))+
  ggtitle("AI-specific technological space (1974-1988)") #
```

We do the same for the 2 other intervals, and combine the three figures again using the multiplot custom function. The resulting figure is saved at `Files_created_with_the_code/figures/Figure_2_ATS_and_AI_core_technologies_3_intervals.jpg`.

# 2. Generating Descriptive Figures

This section details the creation of the paper's descriptive figures. These visualizations illustrate key trends in AI patenting and the evolution of national specialization strategies that motivate our main analysis.

## 2.1. Share of Break-in specialisations (Fig 6 and 7)

Here, we generate the plots showing the share of 'break-in' specializations for each country over time. This metric is central to our paper's narrative and is calculated as the ratio of **coinciding specializations** (specialized in both the general field and its AI application) to the country's total number of **general specializations**. A higher share indicates that a larger portion of a country's established technological strengths is being integrated with AI. We first perform this analysis at the technological field level.


```{r, include=FALSE}
rm(list = ls()[!sapply(ls(), function(x) is.function(get(x)))])
gc()

IPC_RCAs_Top4 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/RCA_4countries_detailed.csv", sep = ";", header = TRUE, dec=",")
IPC_RCAs_Top4$Total_RCA <- as.factor(IPC_RCAs_Top4$Total_RCA)
IPC_RCAs_Top4$Period_sim <- as.numeric(factor(IPC_RCAs_Top4$Period,levels=unique(IPC_RCAs_Top4$Period)))
IPC_RCAs_Top4$techn_field_nr <- as.character(IPC_RCAs_Top4$techn_field_nr)

#replace names:
IPC_RCAs_Top4$ctry_code <- gsub("US", "USA", str_trim(IPC_RCAs_Top4$ctry_code))
IPC_RCAs_Top4$ctry_code <- gsub("CN", "China", str_trim(IPC_RCAs_Top4$ctry_code))
IPC_RCAs_Top4$ctry_code <- gsub("JP", "Japan", str_trim(IPC_RCAs_Top4$ctry_code))
IPC_RCAs_Top4$ctry_code <- gsub("KR", "South Korea", str_trim(IPC_RCAs_Top4$ctry_code))

AI_RCA <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/Specializations_All_periods_IPC.csv", sep = ";", header = TRUE, dec=",")
AI_RCA$Period_sim <- as.numeric(factor(AI_RCA$Period,levels=unique(AI_RCA$Period)))
AI_RCA <- AI_RCA[,c(2,9,13)]
AI_RCA$techn_field_nr <- as.character(AI_RCA$techn_field_nr)
names(AI_RCA) <- c("techn_field_nr", "RCA_AI_Period", "Period_sim")
IPC_RCAs_Top4 <- left_join(IPC_RCAs_Top4, AI_RCA, by = c("techn_field_nr", "Period_sim"))
#fix Total_RCA:
IPC_RCAs_Top4$Total_RCA_2 <- IPC_RCAs_Top4$Round_general + 2*IPC_RCAs_Top4$Round_AI
rm(AI_RCA)

IPC_RCAs_Top4$Coiciding <- ifelse(IPC_RCAs_Top4$Total_RCA_2 ==3,1,0)
IPC_RCAs_Top4$justGeneral <- ifelse(IPC_RCAs_Top4$Total_RCA_2 ==1,1,0)
IPC_RCAs_Top4$OnlyAI <- ifelse(IPC_RCAs_Top4$Total_RCA_2 ==2,1,0)

#now, create a file per country per interval, where I sum over the 3 columns;
IPC_RCAs <- IPC_RCAs_Top4
IPC_RCAs %<>% 
  group_by(ctry_code,Period) %>%   mutate(Share_coinciding = sum(Coiciding)/(sum(Coiciding)+sum(justGeneral))) %>% 
  mutate(Share_OnlyAI = sum(OnlyAI)/(sum(OnlyAI)+sum(Coiciding))) %>% 
  mutate(sum_coinciding = sum(Coiciding)) %>%   mutate(sum_justGeneral = sum(justGeneral)) %>%
  mutate(sum_OnlyAI = sum(OnlyAI)) %>%  ungroup()
```
The data is processed to count the number of 'coinciding', 'general only', and 'AI only' specializations for each country and period. From these counts, the `Share_coinciding` is calculated. The resulting summary table is shown below.
```{r}
SummaryAllData<-distinct(IPC_RCAs, ctry_code, Period, .keep_all = TRUE) 
colnames(SummaryAllData)[1] <- "Country"
head(SummaryAllData)
```
This summarized data is then used to plot the evolution of the break-in share for the four focus countries (Figure 6).

```{r, fig.width=12, fig.height=4}
  ggplot(data=SummaryAllData, aes(x=Period, y=Share_coinciding, group=Country, shape = Country, color=Country)) +
  geom_point(aes(fill = Country), size=8) +   scale_shape_manual(values=c(21, 22, 24, 23)) +
  xlab("Interval") +  ylab("Share of break-in specialisations (%)") +
  theme_classic() +  geom_line(aes(color=Country), linetype = "dashed", size=1.5)+
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A")) +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"))
```

To ensure the robustness of our findings, we repeat the analysis at a more granular level of technological classification: the **4-digit IPC subclass**. This serves as a check to confirm that the observed trends are not an artifact of the broader 35-field aggregation.

```{r, include=FALSE}
IPC_RCAs <- read.csv("Files_created_with_the_code/data/files_code_4-digits_analysis/IPC_RCAs_subclass.csv", sep = ";", header = TRUE, dec=",")
#Select the 4 countries we want
IPC_RCAs_Top4 <- IPC_RCAs[IPC_RCAs$ctry_code == "CN" | IPC_RCAs$ctry_code == "KR"| 
                            IPC_RCAs$ctry_code == "US"|IPC_RCAs$ctry_code == "JP", ]
rm(IPC_RCAs)

#replace names:
IPC_RCAs_Top4$ctry_code <- gsub("US", "USA", str_trim(IPC_RCAs_Top4$ctry_code))
IPC_RCAs_Top4$ctry_code <- gsub("CN", "China", str_trim(IPC_RCAs_Top4$ctry_code))
IPC_RCAs_Top4$ctry_code <- gsub("JP", "Japan", str_trim(IPC_RCAs_Top4$ctry_code))
IPC_RCAs_Top4$ctry_code <- gsub("KR", "South Korea", str_trim(IPC_RCAs_Top4$ctry_code))

IPC_RCAs_Top4$Period_sim <- as.numeric(factor(IPC_RCAs_Top4$Period,levels=unique(IPC_RCAs_Top4$Period)))

#replace NAs by 0:
#replace NAs, so we don't have problems when summing:
IPC_RCAs_Top4[is.na(IPC_RCAs_Top4)] <- 0

#make the numbers binary
IPC_RCAs_Top4$RCA_Gen2 <- ifelse(IPC_RCAs_Top4$RCA_Gen >=1,1,0)
IPC_RCAs_Top4$RCA_AI2 <- ifelse(IPC_RCAs_Top4$RCA_AI >=1,1,0)

#fix Total_RCA:
IPC_RCAs_Top4$Total_RCA_2 <- IPC_RCAs_Top4$RCA_Gen2 + 2*IPC_RCAs_Top4$RCA_AI2

IPC_RCAs_Top4$Coiciding <- ifelse(IPC_RCAs_Top4$Total_RCA_2 ==3,1,0)
IPC_RCAs_Top4$justGeneral <- ifelse(IPC_RCAs_Top4$Total_RCA_2 ==1,1,0)
IPC_RCAs_Top4$OnlyAI <- ifelse(IPC_RCAs_Top4$Total_RCA_2 ==2,1,0)

#now, create a file per country per interval, where I sum over the 3 columns;
IPC_RCAs <- IPC_RCAs_Top4
IPC_RCAs %<>% 
  group_by(ctry_code,Period) %>% 
  mutate(Share_coinciding = sum(Coiciding)/(sum(Coiciding)+sum(justGeneral))) %>% 
  mutate(Share_OnlyAI = sum(OnlyAI)/(sum(OnlyAI)+sum(Coiciding))) %>% 
  mutate(sum_coinciding = sum(Coiciding)) %>% 
  mutate(sum_justGeneral = sum(justGeneral)) %>%
  mutate(sum_OnlyAI = sum(OnlyAI)) %>%
  ungroup()

SummaryAllData4dig<-distinct(IPC_RCAs, ctry_code, Period, .keep_all = TRUE) 
colnames(SummaryAllData4dig)[1] <- "Country"
```
The resulting plot (Figure 7) confirms that the trends observed at the field level are consistent at the more detailed subclass level.
```{r, fig.width=12, fig.height=4}
  ggplot(data=SummaryAllData4dig, aes(x=Period, y=Share_coinciding, group=Country, shape = Country, color=Country)) +
  geom_point(aes(fill = Country), size=8) + 
  scale_shape_manual(values=c(21, 22, 24, 23)) +
  xlab("Interval") +
  ylab("Share of break-in specialisations (%)") +
  theme_classic() +
  geom_line(aes(color=Country), linetype = "dashed", size=1.5)+
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A")) +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A")) 
```

## 2.2. Growth of AI Patents (Fig 1)

This section reproduces Figure 1 from the paper, which illustrates the dramatic growth in AI patenting since the 1970s. We use the raw AI patent data, aggregating the number of unique patent applications per country for each year. A **log-10** scale is used for the y-axis to accommodate the exponential increase in patent counts and allow for a clearer comparison of growth trajectories between the four focus countries, resulting in the figure seen below.

```{r, include=FALSE}
patents_AI_specific <- read.csv("other_files/IPCs_AI.csv", sep = ";", header = TRUE, dec=",")
patents_AI_specific <- patents_AI_specific[,c((1), (3:4))]

patents_AI_specific %<>%  mutate(DistinctOwnerInf = !duplicated(appln_id)) %>%  ungroup()

patents_AI_specific %<>%   group_by(appln_id) %>% 
  mutate(DistinctpatentOffice = !duplicated(patent_office)) %>%  ungroup()

patents_AI_specific %<>%  group_by(appln_id) %>% 
  mutate(DistinctpatentOffice = n_distinct(patent_office, na.rm = T)) %>%  ungroup()

table(patents_AI_specific$DistinctpatentOffice)
test<- patents_AI_specific[1,]

test$patent_office <- gsub("CN", "US", str_trim(test$patent_office))

patents_AI_specific2 <- rbind(patents_AI_specific, test)

patents_AI_specific2 %<>% group_by(appln_id) %>% 
  mutate(DistinctpatentOffice = n_distinct(patent_office, na.rm = T)) %>% ungroup()

patents_AI_specific2[patents_AI_specific2$appln_id == "475222998",]
table(patents_AI_specific2$DistinctpatentOffice)
#thus, there is no patent with inventors from distinct patent offices in our dataset;

patents_AI_specific_simplified <- patents_AI_specific[patents_AI_specific$DistinctOwnerInf == T,]
patents_AI_specific_simplified2 <- patents_AI_specific2[patents_AI_specific2$DistinctOwnerInf == T,]

patents_AI_specific_simplified_4<- patents_AI_specific_simplified[patents_AI_specific_simplified$patent_office == "CN" |
                                                                    patents_AI_specific_simplified$patent_office == "US"|
                                                                    patents_AI_specific_simplified$patent_office == "KR"|
                                                                    patents_AI_specific_simplified$patent_office == "JP", ]

patents_AI_specific_simplified_4$patent_office <- gsub("US", "USA", str_trim(patents_AI_specific_simplified_4$patent_office))
patents_AI_specific_simplified_4$patent_office <- gsub("CN", "China", str_trim(patents_AI_specific_simplified_4$patent_office))
patents_AI_specific_simplified_4$patent_office <- gsub("JP", "Japan", str_trim(patents_AI_specific_simplified_4$patent_office))
patents_AI_specific_simplified_4$patent_office <- gsub("KR", "South Korea", str_trim(patents_AI_specific_simplified_4$patent_office))

table(patents_AI_specific_simplified_4$patent_office)
Data <- as.data.frame(table(patents_AI_specific_simplified_4$patent_office, patents_AI_specific_simplified_4$priority_year))
names(Data) <- c("Country", "Year", "Number_of_AI_patents")

Data$Year <- as.Date(paste(Data$Year, 1, 1, sep = "-")) # beginning of year
Data$Year <- as.Date(paste(Data$Year, 12, 31, sep = "-"))
Data$Year <- as.numeric(format(Data$Year, "%Y"))

Data$Period <- ifelse(Data$Year >= 1974 & Data$Year <= 1988, "First Period (1974-1988)",
                      ifelse(Data$Year > 1988 & Data$Year <= 2003, "Second Period (1989-2003)",
                             ifelse(Data$Year >= 2004 & Data$Year < 2019, "Third Period (2004-2018)", "No period")))
test <- Data[Data$Period != "No period",]  
```

```{r, fig.width=12, fig.height=8}
  ggplot(data=test, aes(x=Year, y=log10(Number_of_AI_patents), group=Country, colour=Country, shape=Country)) +
  geom_line(size=1.2, aes(linetype=Country)) +
  geom_point(size=4) +  xlab("Year") +  ylab("Number of new AI registers [Log10]") + theme_classic() +
  scale_linetype_manual(values=c("twodash", "longdash", "solid", "solid")) +
  scale_shape_manual(values=c(16, 15, 17, 18)) + theme(legend.position="bottom") +
  theme(text = element_text(size = 15)) +  scale_y_continuous(limits=c(0,4)) + 
  geom_vline(data=test, aes(xintercept=c(1988),  colour=Period), linetype="dashed", size=1, color = "grey") +  
  geom_vline(data=test, aes(xintercept=c(2003),  colour=Period), linetype="dashed", size=1, color = "grey") +  
  scale_x_continuous(breaks = c(1974, 1988, 2003, 2018), limits=c(1974, 2018)) + scale_color_brewer(palette="Dark2") + 
  annotate("rect", xmin = 1974, xmax=1988, ymin = 3.6, ymax = 4, alpha = .01, color = "black") +
  annotate("text", x = 1981, y = 3.8, label = c("First Interval \n(1974-1988)"), size=4)+
  annotate("rect", xmin = 1988, xmax=2003, ymin = 3.6, ymax = 4, alpha = .01, color = "black") +
  annotate("text", x = 1996, y = 3.8, label = c("Second Interval \n(1989-2003)"), size=4) +
  annotate("rect", xmin = 2003, xmax=2018, ymin = 3.6, ymax = 4, alpha = .01, color = "black") +
  annotate("text", x = 2011, y = 3.8, label = c("Third Interval \n(2004-2018)"), size=4)
```

# 3. Robustness Checks: Permutation Analysis

To ensure that our findings are statistically robust and not merely the result of random chance, we conduct a permutation analysis. The core idea is to create a "null model" by generating thousands of randomized AI patent datasets. By comparing our actual results to the distribution of results from these random datasets, we can assess the statistical significance of our observations. This section details the creation of these permuted datasets and the subsequent recalculation of specialization metrics.

## 3.1. Permutate the AI dataset

The first step is to generate the randomized, or **permuted**, datasets. For each of the four focus countries and for each time interval, we follow a specific procedure:

1.  Count the number of **actual** AI patents the country has in that interval.
2.  Randomly select the **same number** of patents from that country's **entire pool** of patents (both AI and non-AI) for that interval.
3.  Treat this random sample as the new, 'permuted' AI dataset for that country.

This process is repeated 1,000 times (note: `num_permutations` is set to 10 in this example for faster execution) to create 1,000 counterfactual scenarios where 'AI' patents are just random draws from a country's overall technological portfolio.

We begin by reloading the patent data for the first interval (1974-1988) to establish the pool from which random patents will be drawn.

```{r, include=FALSE}
#### Main Code
# Load required libraries
library(tidyverse)
library(magrittr)
library(tidygraph)
library(ggraph)
library(EconGeo)
library(data.table)
library(netrankr)
library(dplyr)
library(tidyr)
library(ggrepel)
library(scales)
library(patchwork)
library(RColorBrewer)
library(janitor) #also used in the clean_names() function
library(ggforce)
library(stringr)
library(openxlsx)
library(gridExtra) #for grid.arrange
library(readxl) #for reading the xlsx files
library(lmtest) #for LM analysis and robustness econometric test
library(sandwich) 
library(stargazer) #for generating nice econometric tables

rm(list=ls())
#set the working directory to where you saved the R code:
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#4.1.First part: Technological Spaces based on Technological field 
group_by_applnID <- function(data){
  data %>%
    group_by(appln_id) %>%
    mutate(field_weight = 1 / n()) %>%
    ungroup()
}

# This function aggregates the weighted fields at the country-technology field level.
group_by_ctry_and_techn_field <- function(data){
  data %<>%
    group_by(ctry_code, techn_field_nr) %>%
    summarise(n_tech_reg = sum(field_weight)) %>%
    ungroup() %>%
    drop_na() 
}


##1.1. Calculate Specializations for Different Time intervals
###1.1.1 Read and prepare big files
# Reading large files in parts to avoid memory issues
rows_part2 <- 45182803 - 40000000
ipc_all_patents_part2_chunk1_df <- fread("large_files/All_patents_and_IPCs_Part2.csv", 
                                         header = FALSE, nrow = 20000000)
ipc_all_patents_part2_chunk2_df <- fread("large_files/All_patents_and_IPCs_Part2.csv", 
                                         header = FALSE, nrow = 20000000, skip = 20000000)
ipc_all_patents_part2_chunk3_df <- fread("large_files/All_patents_and_IPCs_Part2.csv", 
                                         header = FALSE, nrow = rows_part2, skip = 40000000)

ipc_all_patents_part2_df <- rbind(ipc_all_patents_part2_chunk1_df, 
                                  ipc_all_patents_part2_chunk2_df, ipc_all_patents_part2_chunk3_df)

rm(ipc_all_patents_part2_chunk1_df, ipc_all_patents_part2_chunk2_df, ipc_all_patents_part2_chunk3_df)

setnames(ipc_all_patents_part2_df, c("appln_id", "ctry_code", "techn_field_nr", "weight", "priority_year"))

# Load AI-specific data
ai_patents_df <- fread("other_files/IPCs_AI.csv", sep=";", dec=",", header=TRUE)
ai_patents_df$ctry_code <- "AI_pat"

# Load IPC technology names for labeling
ipc_names_df <- read.csv("other_files/ipc_technology.csv", sep = ";", header = TRUE) %>%
  select(field_nr, sector, field_name) %>% distinct(field_nr, .keep_all = TRUE) %>%
  mutate(techn_field_nr = field_nr) %>% arrange(techn_field_nr)

###1.1.2.Calculate per interval
####1.1.2.1. Interval 1 (1974 - 1988)
# Define Time Interval
start_year <- 1973
end_year   <- 1989

###Filter Data for Interval 1 
ipc_all_patents_first_period_df <- ipc_all_patents_part2_df[priority_year > start_year & priority_year < end_year]
ipc_all_patents_first_period_df <- ipc_all_patents_first_period_df[, .(appln_id, ctry_code, techn_field_nr)]
# Compute weighted values at country-technology field level
region_tech_fields_1_df <- group_by_applnID(ipc_all_patents_first_period_df)
region_tech_fields_1_df <- group_by_ctry_and_techn_field(region_tech_fields_1_df)
```

The resulting count of technological fields per country is:

```{r}
kable(as.data.frame(region_tech_fields_1_df[1:6,]))
```

The following loop executes the permutation logic. For each of the 10 iterations, it samples a new set of random 'AI' patents for the target countries.

```{r, include=FALSE}
ai_patents_df <- fread("other_files/IPCs_AI.csv", sep=";", dec=",", header=TRUE)
ai_patents_period_1_df <- ai_patents_df[ai_patents_df$priority_year > start_year & ai_patents_df$priority_year < end_year,]
# Define target countries
target_countries <- c("CN", "JP", "US", "KR") 
num_permutations <- 10 # Or the desired number
```

```{r}

list_of_permuted_dfs <- vector("list", length = num_permutations)

for (p in 1:num_permutations) {
  if (p %% 100 == 0) print(paste("Permutation number:", p)) # Progress indicator
  
  # This dataframe will hold the permuted AI patents for target countries ONLY for THIS iteration
  permuted_ai_for_target_countries_iter <- data.frame()
  
  for (country in target_countries) {
    # 1. Identify and Count ACTUAL AI patents for the current country from the original AI dataset
    actual_ai_appln_ids_country <- ai_patents_period_1_df %>%
      filter(ctry_code == country) %>%
      distinct(appln_id) %>%
      pull(appln_id)
    
    n_ai_country <- length(actual_ai_appln_ids_country)
    
    if (n_ai_country == 0) {
      # print(paste("No AI patents found for", country, "in original AI data. Skipping for perm", p))
      next # Skip to the next country if no AI patents to replace
    }
    
    # 2. Prepare the pool of ALL patents for the current country from the general dataset
    country_all_patents_pool <- ipc_all_patents_first_period_df %>%
      filter(ctry_code == country) %>%
      distinct(appln_id)
    
    if (nrow(country_all_patents_pool) == 0) {
      # print(paste("No patents in general pool for", country, ". Skipping for perm", p))
      next
    }
    
    # Handle cases where the pool is smaller than the number of AI patents to sample
    # This is unlikely if ipc_all_patents_first_period_df is complete, but good for robustness
    sample_size <- min(n_ai_country, nrow(country_all_patents_pool))
    replace_sampling <- FALSE
    if (n_ai_country > nrow(country_all_patents_pool)) {
      # print(paste("Warning: For country", country, "in perm", p,
      #             "not enough unique patents in pool. Sampling", nrow(country_all_patents_pool),
      #             "instead of", n_ai_country, "OR consider sampling with replacement."))
      # Decide: either sample fewer (as done with min()), or sample with replacement.
      # If sampling with replacement is desired:
      # sample_size <- n_ai_country
      # replace_sampling <- TRUE
      # For now, we sample up to the available pool size without replacement.
      # Or, if strict adherence to n_ai_country is needed and pool is too small WITH replace=FALSE:
      if(nrow(country_all_patents_pool) < n_ai_country && !replace_sampling){
        # print(paste("Strict N_AI needed, but pool too small for", country, "in perm", p, ". Skipping country for this perm."))
        next # Skip this country for this permutation if not enough patents
      }
    }
    
    
    # 3. Randomly select an equivalent number of unique appln_ids from this country's general pool
    random_appln_ids_country <- sample(country_all_patents_pool$appln_id,
                                       size = sample_size, # Use adjusted sample_size
                                       replace = replace_sampling) # Use replace_sampling flag
    
    # 4. Get all rows for these randomly selected patents from the ipc_all_patents_first_period_df
    randomly_selected_patents_df_country <- ipc_all_patents_first_period_df %>%
      filter(appln_id %in% random_appln_ids_country & ctry_code == country)
    
    # 5. Add these randomly selected patents for the current country to the iteration's df
    if (nrow(randomly_selected_patents_df_country) > 0) {
      permuted_ai_for_target_countries_iter <- bind_rows(
        permuted_ai_for_target_countries_iter,
        randomly_selected_patents_df_country
      )
    }
  } # End of country loop
  
  # Add the permutation number to all rows of this iteration's dataframe
  if (nrow(permuted_ai_for_target_countries_iter) > 0) {
    permuted_ai_for_target_countries_iter$permutation_number <- p
  }
  
  # Store the dataframe for this iteration in the list
  list_of_permuted_dfs[[p]] <- permuted_ai_for_target_countries_iter
  
} # End of permutation loop

# Combine all permuted dataframes from the list into one large dataframe
final_permuted_dataset <- bind_rows(list_of_permuted_dfs)
```

The resulting permuted dataset looks like this for the 6 initial and 6 last lines:

```{r}
kable(as.data.frame(final_permuted_dataset[1:6,]))
tail(final_permuted_dataset)
```

This process creates a long-format dataframe where each `permutation_number` represents a complete, unique, randomized AI dataset. A summary for the number of patents for Japan and the US in the first 5 permutations is shown below (only these countries had patents in the first interval; China and South Korea join in the second interval).

```{r}
final_permuted_dataset %>%
  filter(permutation_number <= 5) %>%
  group_by(permutation_number, ctry_code) %>%
  summarise(unique_appln_ids = n_distinct(appln_id), .groups = 'drop') %>%
  print(n=20)
```
A crucial step is to handle the non-target countries. Since our hypothesis is not about them, their actual AI patents are kept constant and are simply replicated across all 1,000 permutations. This ensures that the global context remains stable while only the composition of AI within our focus countries is randomized.

```{r, include=FALSE}
'%notin%' <- Negate('%in%')
not_selected_AI <- ai_patents_period_1_df[ai_patents_period_1_df$ctry_code %notin% target_countries, c("appln_id", "ctry_code")]
list_of_replicated_not_selected_ai <- vector("list", length = num_permutations + 1)
# Loop from 0 to num_permutations
for (i in 0:num_permutations) {
  # Create a copy of the not_selected_AI dataframe for this iteration
  temp_df <- not_selected_AI
  
  # Add the permutation_number column
  temp_df$permutation_number <- i
  
  # Store it in the list (adjust index because list is 1-based, i is 0-based)
  list_of_replicated_not_selected_ai[[i + 1]] <- temp_df
  
  if (i %% 100 == 0) print(paste("Replicating not_selected_AI for permutation_number:", i))
}

# Combine all replicated dataframes from the list into one large dataframe
replicated_not_selected_ai_final <- bind_rows(list_of_replicated_not_selected_ai)
```

The resulting dataset looks like this for the first and last 6 observations:

```{r}
kable(as.data.frame(replicated_not_selected_ai_final[1:6,]))
tail(replicated_not_selected_ai_final)
```

Finally, the permuted data for the target countries is combined with the replicated data for non-target countries. We also add the original, non-permuted AI dataset, labeling it as `permutation_number = 0`. This allows for direct comparison. The entire collection is then joined with the technological field information to prepare for the RTA calculation.

```{r, include=FALSE}
original_selected_AI <- ai_patents_period_1_df[ai_patents_period_1_df$ctry_code %in% target_countries, c("appln_id", "ctry_code")]
table(original_selected_AI$ctry_code)
#table(original_selected_AI$priority_year) #it's correct, so I can exclude this
original_selected_AI$permutation_number <- 0
final_permuted_dataset <- subset(final_permuted_dataset, select = -c(techn_field_nr) )
#rbind it with the original dataset:
final_permuted_dataset <- rbind(original_selected_AI,final_permuted_dataset) #kind of strange and suspicious
table(final_permuted_dataset$permutation_number)
#rbind with the original dataset of NOT SELECTED COUNTRIES
final_permuted_dataset <- rbind(final_permuted_dataset,replicated_not_selected_ai_final)
final_permuted_dataset$ctry_code <- "AI_pat"
# Join technological fields info to AI patents
final_permuted_dataset <- distinct(final_permuted_dataset, appln_id, .keep_all = TRUE)[, c("appln_id", "permutation_number")]
final_permuted_dataset <- left_join(final_permuted_dataset, ipc_all_patents_first_period_df, by = "appln_id")
```

The result is a single dataframe containing the original AI data (permutation 0) and 10 random variations (again, 10 is used for this illustrative example; the number of permutations was set to 1000 in the files that are saved in the folder `Files_created_with_the_code/data/files_code_Fields_analysis/robustness/`). A summary of this file for this illustrative example looks like this:

```{r}
table(final_permuted_dataset$permutation_number)
```

## 3.2. Calculate AI-specific specialisations

With the 11 datasets (1 actual + 10 permuted) assembled, we now repeat the exact same **Revealed Technological Advantage (RTA)** calculation performed in Section 1.1. This is done for each permutation, allowing us to generate a distribution of RTA scores for each technological field under the null hypothesis (i.e., when 'AI' is random).

```{r, include=FALSE}
group_by_applnID_perm <- function(data){
  # No change needed here if field_weight is per appln_id irrespective of permutation
  # However, if an appln_id could theoretically appear in multiple permutations
  # (which it shouldn't if you sampled unique appln_ids per permutation for target countries,
  # and replicated non-target countries correctly), this function is fine as is.
  # The grouping is on appln_id itself.
  data %>%
    group_by(appln_id, permutation_number) %>% # Add permutation_number if an appln_id might exist across permutations
    # If appln_id is unique within each permutation, then just appln_id is fine
    mutate(field_weight = 1 / n()) %>%
    ungroup()
}

group_by_ctry_techn_field_perm <- function(data){
  data %>% # Removed %<>% as it's better to assign explicitly in the main flow
    group_by(permutation_number, ctry_code, techn_field_nr) %>% # Add permutation_number
    summarise(n_tech_reg = sum(field_weight), .groups = 'drop') %>% # Use .groups = 'drop'
    # ungroup() %>% # .groups = 'drop' handles this
    drop_na()
}
region_tech_fields_perm_df <- group_by_applnID_perm(final_permuted_dataset)
region_tech_fields_perm_df <- group_by_ctry_techn_field_perm(region_tech_fields_perm_df)
```
The code below iterates through each `permutation_number`, calculates the AI-specific RTAs for that dataset, and stores the results.
```{r}
list_of_rca_dfs <- region_tech_fields_perm_df %>%
  group_by(permutation_number) %>%
  group_split() %>% # This splits the df into a list of dfs, one for each permutation
  purrr::map(~{
    current_permutation_number <- unique(.x$permutation_number)
    print(paste("Processing RCA for permutation_number:", current_permutation_number))
    
    # Matrix creation for the current permutation's data
    mat_reg_tech_perm_AI <- .x %>%
      select(-permutation_number) %>% # Temporarily remove for pivot if it causes issues
      arrange(techn_field_nr, ctry_code) %>%
      pivot_wider(names_from = techn_field_nr,
                  values_from = n_tech_reg,
                  values_fill = 0) # Changed from list(n_tech_reg = 0) for simplicity
    
    # Check if ctry_code column exists and is not empty
    if (!"ctry_code" %in% names(mat_reg_tech_perm_AI) || nrow(mat_reg_tech_perm_AI) == 0 || all(is.na(mat_reg_tech_perm_AI$ctry_code))) {
      print(paste("Skipping permutation", current_permutation_number, "due to missing ctry_code or empty data after pivot."))
      return(NULL) # Return NULL or an empty tibble
    }
    
    # Check for duplicate ctry_codes which would prevent rownames_to_column
    if (any(duplicated(mat_reg_tech_perm_AI$ctry_code))) {
      print(paste("Warning: Duplicate ctry_code found for permutation", current_permutation_number, ". Aggregating or handling needed."))
      return(tibble(permutation_number = current_permutation_number, error="duplicate ctry_code"))
    }
    
    
    mat_reg_tech_perm_AI <- mat_reg_tech_perm_AI %>%
      remove_rownames() %>%
      column_to_rownames(var = "ctry_code") %>%
      as.matrix() %>%  round()# No rounding here, location_quotient might prefer raw numbers
    
    # RCA calculation
    # Ensure matrix is suitable (e.g., no NA/NaN/Inf that location_quotient can't handle)
    if (nrow(mat_reg_tech_perm_AI) == 0 || ncol(mat_reg_tech_perm_AI) == 0) {
      print(paste("Skipping RCA for permutation", current_permutation_number, "due to empty matrix."))
      return(NULL)
    }
    
    # Ensure there are at least two columns for location_quotient (ctry_code was one)
    if (ncol(mat_reg_tech_perm_AI) < 1) { # If only ctry_code was present and now it's rownames
      print(paste("Skipping RCA for permutation", current_permutation_number, "due to insufficient columns in matrix."))
      return(NULL)
    }
    
    
    # Check for all zero rows/columns if location_quotient is sensitive
    # For example, if a row sum is 0, RCA might be NaN or Inf.
    # The location_quotient function might handle this, or you might need pre-filtering.
    
    rca_results_perm <- tryCatch({
      mat_reg_tech_perm_AI %>%
        location_quotient(binary = FALSE) %>% 
        as.data.frame() %>%
        rownames_to_column("ctry_code") %>%
        as_tibble() %>%
        gather(key = "techn_field_nr", value = "RCA", -ctry_code) %>%
        arrange(ctry_code, techn_field_nr) %>%
        mutate(permutation_number = current_permutation_number) # Add back permutation number
    }, error = function(e) {
      print(paste("Error in location_quotient for permutation", current_permutation_number, ":", e$message))
      return(tibble(permutation_number = current_permutation_number, ctry_code=NA, techn_field_nr=NA, RCA=NA, error_message = e$message)) # Return an empty or error-marked tibble
    })
    
    return(rca_results_perm)
  })

# Combine the list of RCA dataframes into one final dataframe
final_rca_all_permutations_df <- bind_rows(list_of_rca_dfs)
```
The final output is a comprehensive dataframe containing the calculated RTA scores for every country, technology, and permutation, which looks like this for the first and last 6 rows:

```{r}
kable(as.data.frame(final_rca_all_permutations_df[1:6,]))
tail(final_rca_all_permutations_df)
```

This process is repeated for all three time intervals, and the resulting dataframes are saved at `Files_created_with_the_code/data/files_code_Fields_analysis/robustness/`. These files form the basis for the statistical tests in our econometric analysis. It is also worth noting that these permutations are repeated for several distinct interval lengths, i.e., 15-years, 10-years, 5-years, and 1-year.

# 4. Regression Analysis

This final section presents the econometric analysis designed to formally test our paper's hypotheses. Using the data prepared in the previous steps, we construct a panel dataset covering our four focus countries across nine 5-year intervals (from 1974-1978 to 2014-2018). We then run a series of regression models to investigate the factors influencing the emergence and persistence of different types of technological specializations.

The following code block handles the final data preparation, loading the pre-calculated metrics for relative density and specializations, and merging them into a single dataframe ready for regression.


```{r, include=FALSE}
options(scipen = 999) #deactivate scientific notation
rm(list=ls())
Rel_density <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness/Rel_density_5y.csv", 
                        sep = ";", header = TRUE, dec=",")
target_countries <- c("CN", "JP", "US", "KR") 
Rel_density <-  Rel_density[Rel_density$ctry_code %in% target_countries,] 

Periods_5y <- c("1974-1978","1979-1983","1984-1988","1989-1993","1994-1998","1999-2003","2004-2008","2009-2013","2014-2018")

Per_1 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_1974-1978.csv", 
                 sep = ";", header = TRUE, dec=",")
Per_2 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_1979-1983.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_3 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_1984-1988.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_4 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_1989-1993.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_5 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_1994-1998.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_6 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_1999-2003.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_7 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_2004-2008.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_8 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_2009-2013.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_9 <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness_automated/Data_RCA_techn_field_5_years_2014-2018.csv", 
                  sep = ";", header = TRUE, dec=",")
Per_all <- rbind(Per_1, Per_2, Per_3, Per_4, Per_5, Per_6, Per_7, Per_8, Per_9)
rm(Per_1, Per_2, Per_3, Per_4, Per_5, Per_6, Per_7, Per_8, Per_9)
unique(Per_all$Period)

Per_all[is.na(Per_all)] <- 0
Per_all <- Per_all[Per_all$ctry_code %in% target_countries,] 
Per_all$Round_general <- ifelse(Per_all$RCA_Gen < 1, 0, 1)
Per_all$Round_AI <- ifelse(Per_all$RCA_AI < 1, 0, 1)
Per_all$Total_RCA_2 <- Per_all$Round_general + 2*Per_all$Round_AI

Specializations_data <- as.data.frame(table(Per_all$Total_RCA_2, Per_all$ctry_code, Per_all$Period))
Specializations_data$Var1 <- gsub("0", "No specialization", str_trim(Specializations_data$Var1))
Specializations_data$Var1 <- gsub("1", "General specialization", str_trim(Specializations_data$Var1))
Specializations_data$Var1 <- gsub("2", "AI-specific specialization", str_trim(Specializations_data$Var1))
Specializations_data$Var1 <- gsub("3", "Coinciding specialization", str_trim(Specializations_data$Var1))

names(Specializations_data) <- c("Var1","ctry_code", "Period", "N_spec")
# Pivot the data
Specializations_wide <- Specializations_data %>%
  pivot_wider(names_from = Var1,
              values_from = N_spec)

Rel_density <- left_join(Rel_density, Specializations_wide, by = c("ctry_code", "Period"))

actual_shares_5_years <- read.csv("Files_created_with_the_code/data/files_code_Fields_analysis/robustness/actual_shares_5_years.csv", 
                                   sep = ";", header = TRUE, dec=",")
Rel_density <- left_join(Rel_density, actual_shares_5_years, by = c("ctry_code", "Period"))
Rel_density %<>%  clean_names()
Rel_density$total_general_specializations <- Rel_density$general_specialization + Rel_density$coinciding_specialization
Rel_density$double_check <- Rel_density$total_general_specializations + Rel_density$no_specialization + Rel_density$ai_specific_specialization
Rel_density$total_specializations <- Rel_density$total_general_specializations + Rel_density$ai_specific_specialization #which should be the same as
#total_general_specializations - no_specialization; pay attention that this includes AI
Rel_density$total_AI_specializations <- Rel_density$coinciding_specialization + Rel_density$ai_specific_specialization

regression_data <- Rel_density %>%
  mutate(Share_Coinciding =  coinciding_specialization/(total_general_specializations),
         # Ensure Period is a factor for the regression
         Period = factor(period, levels = c("1974-1978","1979-1983","1984-1988","1989-1993","1994-1998","1999-2003",
                                            "2004-2008","2009-2013","2014-2018"), ordered = FALSE), # Not ordered for lm factor
         ctry_code = factor(ctry_code))
rm(Specializations_wide, Specializations_data)

regression_data <- regression_data %>%
  mutate(ctry_code = case_when(
    ctry_code == "JP" ~ "Japan",
    ctry_code == "US" ~ "US",
    ctry_code == "KR" ~ "South Korea",
    ctry_code == "CN" ~ "China",
    TRUE ~ ctry_code # This keeps any other codes as they are
  ))

#redefine names:
regression_data_renamed <- regression_data %>%
  rename(
    "Country’s techn. rel. dens." = rel_density,
    "Interval " = Period,
    "No. of ‘general’ spec." = total_general_specializations,
    "Share of ‘break-in’ spec." = Share_Coinciding,
    "No. of ‘AI-specific’ spec." = total_AI_specializations,
    "No. of sustained ‘break-in’ spec." = actual_persistent_coinciding,
    "No. of sustained ‘general’ spec." = actual_persistent_general_all,
    "No. of sustained ‘AI-specific’ spec." = actual_n_persistent_round_ai,
    "No. of ‘break-in’ spec." = coinciding_specialization,
    "Country "= ctry_code
  )

desired_order_model1_to_3 <- c("`Country’s techn. rel. dens.`", "`Interval `", "`Country `", "`No. of ‘general’ spec.`",
                   "`No. of ‘AI-specific’ spec.`")

# Model 1: Predicting Share_Coinciding with Period fixed effects
model1 <- lm(`Share of ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `, 
             data = regression_data_renamed) 
summary(model1)

model2 <- lm(`Share of ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `No. of ‘general’ spec.` + `Interval ` + 
               `No. of ‘AI-specific’ spec.`, 
             data = regression_data_renamed) 
summary(model2)

# Model 3: A simpler model if degrees of freedom are a major concern for Model 2
model3 <- lm(`Share of ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `No. of ‘general’ spec.` + `Interval ` + `Country ` + 
               `No. of ‘AI-specific’ spec.`,
             data = regression_data_renamed)
summary(model3)
```

The final dataset for regression looks like this:
```{r}
head(regression_data_renamed) %>%
  knitr::kable(
    caption = "Preview of the Final Regression Dataset",
    booktabs = TRUE # A style option for prettier tables
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE
  )
```

## 4.1. Main models (i.e., the ones from the paper) 

Our main analysis consists of two sets of Ordinary Least Squares (OLS) models. The first set (Table 3) examines the factors that influence a country's **share of 'break-in' specializations**. The second set (Table 4) investigates the determinants of the **persistence of specializations over time**.

The first three models test the effect of technological relatedness density on the share of break-in specializations. Model 1 provides a baseline, Model 2 adds control variables for the number of general and AI-specific specializations, and Model 3 includes country fixed effects. The estimations for the three first models are:

```{r, results='asis', echo=FALSE, warning=FALSE, message=FALSE}
stargazer(model1, model2, model3, order = desired_order_model1_to_3, title="Effects on the share of break-ins - OLS regression", ci.level=0.95, single.row=TRUE, ci=F, type="html") 
```

The next set of models shifts the focus to persistence. The dependent variable is now the count of specializations of a certain type that are sustained from the previous period. These models help us understand what factors contribute to the durability of a country's technological advantages. The estimations for these three models are:

```{r, include=FALSE}
newmodel4 <- lm(`No. of sustained ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of sustained ‘general’ spec.` + 
                  `No. of sustained ‘AI-specific’ spec.`, data = regression_data_renamed)
summary(newmodel4)

#Model 5 - persistent break-ins
newmodel5 <- lm(`No. of sustained ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of ‘general’ spec.` + `No. of ‘break-in’ spec.` +
                  `No. of ‘AI-specific’ spec.`, 
                data = regression_data_renamed)
summary(newmodel5)

#Model 6 persistent AI-specific
newmodel6 <- lm(`No. of sustained ‘AI-specific’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of ‘general’ spec.` +
                  `No. of sustained ‘break-in’ spec.`+`No. of ‘break-in’ spec.`+`No. of ‘AI-specific’ spec.`, 
                data = regression_data_renamed)
summary(newmodel6)

desired_order_model4_to_6 <- c("`Country’s techn. rel. dens.`", "`Interval `", 
                               "`No. of ‘break-in’ spec.`", "`No. of ‘general’ spec.`", "`No. of ‘AI-specific’ spec.`",
                               "`No. of sustained ‘break-in’ spec.`", "`No. of sustained ‘general’ spec.`", 
                               "`No. of sustained ‘AI-specific’ spec.`")
```

```{r, results='asis', echo=FALSE, warning=FALSE, message=FALSE}
stargazer(newmodel5, newmodel4,newmodel6, order = desired_order_model4_to_6, title="Effects on persisting specialisations - OLS regression", ci.level=0.95, single.row=TRUE, ci=F, type="html") 
```

## 4.2. Extensions (i.e., models used for additional robustness, not included in the paper) 

To ensure the robustness of our results, we re-estimate our models using specifications better suited to the nature of our dependent variables.

- For the models predicting the **share** of break-ins (a value between 0 and 1), we use **Beta regression**.

- For the models predicting the **count** of persistent specializations, we use **Poisson** and **Negative Binomial regressions**, which are designed for count data.

The Beta regression results for the share of break-ins are presented below.
```{r, include=FALSE}
library(betareg)
library(MASS)
library(statmod)
regression_data_renamed$`Share of ‘break-in’ spec.` <- ifelse(regression_data_renamed$`Share of ‘break-in’ spec.` ==0,0.01, regression_data_renamed$`Share of ‘break-in’ spec.`)
table(regression_data_renamed$`Share of ‘break-in’ spec.`)

# Model 1 - The simplest model
model1_beta <- betareg(`Share of ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `, 
                       data = regression_data_renamed) 
summary(model1_beta)

# Model 2 - adding number of general specialisations
model2_beta <- betareg(`Share of ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `No. of ‘general’ spec.` + `Interval ` + 
                         `No. of ‘AI-specific’ spec.`, 
                       data = regression_data_renamed) 
summary(model2_beta)

# Model 3 - controlling for country
model3_beta <- betareg(`Share of ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `No. of ‘general’ spec.` + `Interval ` + `Country ` + 
                         `No. of ‘AI-specific’ spec.`,
                       data = regression_data_renamed)
summary(model3_beta)

```

```{r, results='asis', echo=FALSE, warning=FALSE, message=FALSE}
stargazer(model1_beta, model2_beta, model3_beta, order = desired_order_model1_to_3, title="Effects on the share of break-ins - Beta regression", ci.level=0.95, single.row=TRUE, ci=F, type="html") 
```
Compared to the OLS-based Table 3, the estimations from the Beta regression are highly consistent. The main differences are minor shifts in significance levels for some time intervals and the US country dummy. Importantly, the main variable of interest, "Country's techn. rel. dens.", retains its significance and sign, confirming the robustness of our primary finding.

Next, we test the robustness of the persistence models. The results from the **Poisson regression** are shown first, followed by the **Negative Binomial models**, which can be more appropriate if the count data is over-dispersed.
  
```{r, include=FALSE}
model4_poisson <- glm(`No. of sustained ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of sustained ‘general’ spec.` + 
                        `No. of sustained ‘AI-specific’ spec.`,
                      family = poisson(link = "log"),
                      data = regression_data_renamed)
summary(model4_poisson)

#Model 5 - persistent break-ins
model5_poisson <- glm(`No. of sustained ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of ‘general’ spec.` + `No. of ‘break-in’ spec.` +
                        `No. of ‘AI-specific’ spec.`, 
                      family = poisson(link = "log"),
                      data = regression_data_renamed)
summary(model5_poisson)

#Model 6 persistent AI-specific
model6_poisson <- glm(`No. of sustained ‘AI-specific’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of ‘general’ spec.` +
                        `No. of sustained ‘break-in’ spec.`+`No. of ‘break-in’ spec.`+`No. of ‘AI-specific’ spec.`, 
                      family = poisson(link = "log"),
                      data = regression_data_renamed)
summary(model6_poisson)
```

```{r, results='asis', echo=FALSE, warning=FALSE, message=FALSE}
stargazer(model5_poisson, model4_poisson,  model6_poisson, order = desired_order_model4_to_6, title="Effects on persisting specialisations - Poisson regression", ci.level=0.95, single.row=TRUE, ci=F, type="html") 
```

For the **Negative binomial**, the results look like this:
```{r, include=FALSE}
model4_negbin <- glm.nb(`No. of sustained ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of sustained ‘general’ spec.` + 
                          `No. of sustained ‘AI-specific’ spec.`,
                        data = regression_data_renamed)
summary(model4_negbin)

#Model 5 - persistent break-ins
model5_negbin <- glm.nb(`No. of sustained ‘break-in’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of ‘general’ spec.` + `No. of ‘break-in’ spec.` +
                          `No. of ‘AI-specific’ spec.`, 
                        data = regression_data_renamed)
summary(model5_negbin)

#Model 6 persistent AI-specific
model6_negbin <- glm.nb(`No. of sustained ‘AI-specific’ spec.` ~ `Country’s techn. rel. dens.` + `Interval `+ `No. of ‘general’ spec.` +
                          `No. of sustained ‘break-in’ spec.`+`No. of ‘break-in’ spec.`+`No. of ‘AI-specific’ spec.`, 
                        data = regression_data_renamed)
summary(model6_negbin)
```

```{r, results='asis', echo=FALSE, warning=FALSE, message=FALSE}
stargazer(model5_negbin, model4_negbin,  model6_negbin, order = desired_order_model4_to_6, title="Effects on persisting specialisations - Negative binomial regression", ci.level=0.95, single.row=TRUE, ci=F, type="html") 
```
The results from both the Poisson and Negative Binomial regressions are nearly identical to each other and largely consistent with the OLS models. While some variables with marginal significance in the OLS models lose significance here, the main relationships of interest hold. This confirms that our findings regarding the persistence of specializations are not sensitive to the choice of a linear versus a count data model specification.