Nuclear Threat Projection with Neural Network Time Series Forecasting

[This article was first published on DataGeeek, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Unfortunately, we have been through tough times recently as going on Russian invasion in Ukraine. As Putin stacked to the corner via sanctions and lost in the field, he has been getting to be more dangerous. He has even threatened to use nuclear weapons if necessary.

Because of the nuclear danger we’ve just mentioned above, I have wondered about the number of nuclear weapons the countries have nowadays and especially the change from the end of the cold war.

The dataset we are going to work on is from Our World in Data. In terms of the easiness of writing, we will change the names of some of the countries.

library(tidyverse)
library(fable)
library(sysfonts)
library(showtext)
library(plotly)
library(readxl)
library(glue)

df <- 
  read_excel("nuclear.xlsx") %>% 
  mutate(country=case_when(
    country == "United Kingdom" ~ "UK",
    country == "United States" ~ "USA",
    country == "North Korea" ~ "PRK",
    country == country ~ country))

We will subset the dataset to dates that we want to examine are 1991, the end of the cold war, and the present day.

#Tidying the dataset for the plot
df_tidy <- 
  df %>% 
  filter(year==1991 | year==2022) %>%
  pivot_wider(names_from = year, values_from = stockpile) %>% 
  mutate(
    color= case_when (
    `1991` - `2022` > 0 ~ "darkgrey",
    `1991` - `2022` < 0 ~ "red"
  )) %>% 
  na.omit() %>% 
  pivot_longer(c(`1991`, `2022`), 
               names_to = "year", 
               values_to = "stockpile") %>%
  mutate(stockpile_pct=percent_rank(stockpile))

Now, we can draw the plot which shows the changes that countries in terms of the number of nuclear warheads at interested years. We will do that with ggplotly as dynamic plotting, so when we hoover the cursor on the plot and country names, we will be able to see the corresponding numbers. We will indicate red the ones on rising and dark grey the ones on the descent.

#Comparing the change of the numbers of nuclear wearheads 

#load fonts(google)
font_add_google("Roboto Slab")
showtext_auto()

p <- 
  df_tidy %>% 
  ggplot(aes(x = year, 
             y = stockpile_pct, 
             group = country, 
             color = color, 
             fill = color,
             text=stockpile)) +
  geom_line(size= 1.5)+
  geom_point(size = 5,stroke = 1) +
  geom_text(
    data= filter(df_tidy,year == "2022"), 
    aes(label = country),
    hjust = c(0.1),
    nudge_x = c(0.18),
    family = "Roboto Slab",
    size = 6
  ) +
  coord_cartesian(clip = "off")+
  scale_y_continuous(limits = c(0, 1), 
                     breaks = scales::pretty_breaks(n = 5))+
  scale_color_identity()+
  scale_fill_identity()+
  scale_size_identity()+
  theme_minimal() + 
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.y = element_blank(),
    axis.text.x = element_text(family = "Roboto Slab", size=20),
    axis.title = element_blank(),
    panel.background = element_rect(fill = "lightblue",color = NA)
  )

#setting font family for ggplotly 
font <- list(family= "Roboto Slab")

ggplotly(p, tooltip = "text") %>% 
  hide_legend() %>% 
  layout(font=font)

It’s interesting that Pakistan and India have entered the arms race since 1991 from zero point. Also, we can clearly see that while Russia and USA have decreased the numbers of warheads, China has increased the stockpiles.

Another important issue is how the stockpiles of each nuclear power country will change over the next decade. Hence, we will make a function that forecasts the next 10-years and shows the results on a plot for each country.

The algorithm we will use is the neural network autoregression (NNAR), and we prefer the function NNETAR from the fable package to implement it. We use the ggplotly function as we used before, again.

#The function that predicts 10-years ahead for each country.
fun_nuc <- 
  function(...){
  
    #A function that automatically quotes the arguments for easy writing.
    #(source: https://adv-r.hadley.nz/quasiquotation.html#quasi-motivation)
    cement <- function(...) {
    args <- ensyms(...)
    paste(purrr::map(args, as_string), collapse = " ")
  }
  
  var_loc <- cement(...) #country value in quotation
  
  #Modeling
  set.seed(1234)
  
  fit <- 
    df %>% 
    filter(country==var_loc) %>% 
    as_tsibble(index = "year") %>%
    select(year, stockpile) %>% 
    model(NNETAR(stockpile ~ year)) 
  
  #Forecasting with NNAR
  pred <- 
    fit %>% 
    forecast(h=10)
  
  #Plot variable
  p_pred <- 
    df[df$country==var_loc,] %>% 
    ggplot(aes(x=year,y=stockpile))+
    geom_line(size=1)+
    geom_point(data=pred, 
               mapping = aes(y=round(.mean)),
               color=df_tidy[["color"]][df_tidy$country == var_loc][1])+
    scale_x_continuous(breaks = scales::pretty_breaks(n=7))+
    scale_y_continuous(breaks = scales::pretty_breaks())+
    xlab("")+
    ylab("")+
    labs(title = var_loc)+
    theme_minimal()+
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(size=15,
                                    colour = df_tidy[["color"]][df_tidy$country== var_loc][1], 
                                    hjust = 0.5))
  
  #Interactive plotting
  ggplotly(p_pred) %>% 
    style(
      text=glue("stockpile: {round(pred$.mean)}\nyear: {pred$year}"),
      traces = 2 #for the predicted time series
    ) %>% 
    layout(font = font)
  
  
}


fun_nuc(China)
fun_nuc(USA)
fun_nuc(Russia)

If we can hover on the line and points, we could see the corresponding year and the numbers of stockpiles thanks to the interactive plotting. It looks like the nuclear danger will move to Asia-Pacific when the crisis is over on the East-Europa; if we would survive.

To leave a comment for the author, please follow the link and comment on their blog: DataGeeek.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)