library(forecast)
library(lmtest)
library(ggplot2)
library(tidyr)

# Reading data
data <- read.csv("/Users/sascha/Nextcloud/R-Scripts/Jackson_Reply_2026/ngram_word2vec_revision_Nov14.csv")

# Select year range
data <- data[data$Year %in% 1961:2014,]

Preparing variables

## Linear interpolation for NA-values in $Patents and $Trademarks ----
### ($Trademarks has no NAs in year range 1961-2014)

### Linear interpolation 
interpolated <- approx(data$Year, data$Patents, xout = data$Year, method = "linear")
### Check interpolation
data.frame(iYear = interpolated$x, oYear = data$Year, iPatents = interpolated$y, oPatents = data$Patents)
##    iYear oYear iPatents oPatents
## 1   1961  1961  79721.0    79721
## 2   1962  1962  83396.0    83396
## 3   1963  1963  84632.5       NA
## 4   1964  1964  85869.0    85869
## 5   1965  1965  87592.0    87592
## 6   1966  1966  94629.0    94629
## 7   1967  1967  88525.0    88525
## 8   1968  1968  88164.0    88164
## 9   1969  1969  93471.0    93471
## 10  1970  1970 101415.0   101415
## 11  1971  1971 103175.0   103175
## 12  1972  1972 104729.0   104729
## 13  1973  1973  99298.0    99298
## 14  1974  1974 104079.0   104079
## 15  1975  1975 102538.0   102538
## 16  1976  1976 101014.0   101014
## 17  1977  1977 102344.0   102344
## 18  1978  1978 100931.0   100931
## 19  1979  1979 100916.0   100916
## 20  1980  1980 100494.0   100494
## 21  1981  1981 104329.0   104329
## 22  1982  1982 106413.0   106413
## 23  1983  1983 109625.0   109625
## 24  1984  1984 103703.0   103703
## 25  1985  1985 111284.0   111284
## 26  1986  1986 115235.0   115235
## 27  1987  1987 120916.0   120916
## 28  1988  1988 131837.0   131837
## 29  1989  1989 143836.0   143836
## 30  1990  1990 158707.0   158707
## 31  1991  1991 171163.0   171163
## 32  1992  1992 172115.0   172115
## 33  1993  1993 183347.0   183347
## 34  1994  1994 184196.0   184196
## 35  1995  1995 202755.0   202755
## 36  1996  1996 228142.0   228142
## 37  1997  1997 211946.0   211946
## 38  1998  1998 220496.0   220496
## 39  1999  1999 236979.0   236979
## 40  2000  2000 265763.0   265763
## 41  2001  2001 295895.0   295895
## 42  2002  2002 326471.0   326471
## 43  2003  2003 334445.0   334445
## 44  2004  2004 342441.0   342441
## 45  2005  2005 356943.0   356943
## 46  2006  2006 390733.0   390733
## 47  2007  2007 425966.0   425966
## 48  2008  2008 456154.0   456154
## 49  2009  2009 456321.0   456321
## 50  2010  2010 456106.0   456106
## 51  2011  2011 490226.0   490226
## 52  2012  2012 503582.0   503582
## 53  2013  2013 542815.0   542815
## 54  2014  2014 571612.0   571612
### Interpolated Patents to data
data$iPatents <- interpolated$y
### z-scores for interpolated Patents
data$ziPatents <- scale(data$iPatents)
### z-scores for Trademarks
data$zTrademarks <- scale(data$Trademarks)

# Just to make sure: Sort by Year
data <- data[order(data$Year),]

Check variables

data[,c("Year", "ziPatents", "zTrademarks")]
##     Year   ziPatents zTrademarks
## 162 1961 -0.87929791 -0.99926584
## 163 1962 -0.85376615 -0.98986275
## 164 1963 -0.84517567 -0.97720004
## 165 1964 -0.83658519 -0.98414198
## 166 1965 -0.82461479 -0.97302923
## 167 1966 -0.77572582 -0.96527003
## 168 1967 -0.81813285 -0.95316155
## 169 1968 -0.82064087 -0.95007102
## 170 1969 -0.78377093 -0.93729559
## 171 1970 -0.72858065 -0.90858844
## 172 1971 -0.71635320 -0.90020925
## 173 1972 -0.70555691 -0.90520670
## 174 1973 -0.74328833 -0.88476601
## 175 1974 -0.71007273 -0.88993254
## 176 1975 -0.72077870 -0.90035016
## 177 1976 -0.73136657 -0.88849531
## 178 1977 -0.72212650 -0.85810668
## 179 1978 -0.73194320 -0.81250026
## 180 1979 -0.73204741 -0.79397588
## 181 1980 -0.73497922 -0.77954715
## 182 1981 -0.70833588 -0.77329094
## 183 1982 -0.69385746 -0.72000987
## 184 1983 -0.67154236 -0.61446223
## 185 1984 -0.71268496 -0.69149059
## 186 1985 -0.66001660 -0.62521802
## 187 1986 -0.63256736 -0.60570730
## 188 1987 -0.59309908 -0.56272172
## 189 1988 -0.51722635 -0.54356796
## 190 1989 -0.43386430 -0.47731418
## 191 1990 -0.33054927 -0.32648890
## 192 1991 -0.24401224 -0.01701339
## 193 1992 -0.23739830 -0.05484182
## 194 1993 -0.15936492 -0.01240108
## 195 1994 -0.15346657  0.19972750
## 196 1995 -0.02452946  0.29963892
## 197 1996  0.15184460  0.56073687
## 198 1997  0.03932424  0.78299186
## 199 1998  0.09872465  0.99059267
## 200 1999  0.21323892  1.10332649
## 201 2000  0.41321336  1.23624737
## 202 2001  0.62255291  1.53405591
## 203 2002  0.83497712  0.81866914
## 204 2003  0.89037582  0.78418486
## 205 2004  0.94592736  0.86241561
## 206 2005  1.04667880  1.12018818
## 207 2006  1.28143200  1.27146436
## 208 2007  1.52621033  1.39423066
## 209 2008  1.73593894  1.64363344
## 210 2009  1.73709916  1.54914220
## 211 2010  1.73560546  1.29339866
## 212 2011  1.97265131  1.43412571
## 213 2012  2.06544101  1.66166935
## 214 2013  2.33800900  1.73298634
## 215 2014  2.53807375  2.00207833
data %>% pivot_longer(cols = c(ziPatents, zTrademarks), # Long format for plotting only
                      names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = Year, y = value, col = variable)) +
  geom_line(alpha = .5) +
  geom_point() +
  theme_minimal()

ARIMAX models

Note: The following defaults apply:

## Outcome: ziPatents; xreg: zTrademarks
Patents_arima <- auto.arima(data$ziPatents, xreg = data$zTrademarks, stepwise = F, seasonal = F)
Patents_arima
## Series: data$ziPatents 
## Regression with ARIMA(0,0,4) errors 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4    xreg
##       1.1265  1.1482  0.8944  0.6593  0.6925
## s.e.  0.1910  0.2686  0.2175  0.1461  0.1377
## 
## sigma^2 = 0.02171:  log likelihood = 27.72
## AIC=-43.44   AICc=-41.65   BIC=-31.5
## Outcome: zTrademarks; xreg: ziPatents
Trademarks_arima <- auto.arima(data$zTrademarks, xreg = data$ziPatents, stepwise = F, seasonal = F)
Trademarks_arima
## Series: data$zTrademarks 
## Regression with ARIMA(0,1,0) errors 
## 
## Coefficients:
##         xreg
##       0.5457
## s.e.  0.1928
## 
## sigma^2 = 0.02344:  log likelihood = 24.77
## AIC=-45.54   AICc=-45.3   BIC=-41.59