Forecasting nuclear proliferation/Simulating nuclear proliferation
Most of the statistical computations summarized in the Wikiversity article on "Forecasting nuclear proliferation" can be replicated by "knitting" the R Markdown code in this vignette in Rstudio.
Process
[edit | edit source]
- Open an instance of RStudio. If you don't already have it installed, you can go to "RStudio Cloud" and click "get started for free" (at least as of 2020-10-24).[1] Or you can download and install free and open-source versions of it and R on your local computer. R is available from The Comprehensive R Archive Network. For RStudio, go to its website > Products (top center) > RStudio > "RStudio Desktop" > "Open Source Edition" > "Download RStudio Desktop".
- Start RStudio. Then File > "New File" > "R Markdown..." > Title: "Modeling and simulating nuclear proliferation" > "Default Output Format: PDF" > OK. [NOTE: The default "Default Output Format" is HTML. When this vignette was developed, the figure references would not work with HTML but would with PDF.]
- Replace the default code in this new RMarkdown vignette on "Modeling and simulating nuclear proliferation" with the text from the section below entitled, 'RMarkdown vignette on "Modeling and simulating nuclear proliferation"'.
- Save as, e.g., "FcstgNucProlif.Rmd" (with a *.Rmd extension).
- File > "New File" > "Text File".
- Copy the text from the section below entitled "nuc-references.bib" into this new empty text file.
- Save as "nuc-references.bib".
- Click the "Knit" icon; see the companion image. Or read the text and run the code chunks one at a time manually. The latter makes it relatively easy to look at intermediate computations carefully and experiment with changing things in different ways.
Development of this vignette
[edit | edit source]This vignette evolved over time and was used to produce most of the figures in the main article on "Forecasting nuclear proliferation" as well as slides for companion presentations. Some were written as svg files. Others were png. These image files were developed in in different formats at different times for different purposes.
RMarkdown vignette on "Modeling and simulating nuclear proliferation"
[edit | edit source]---
title: "Forecasting Nuclear Proliferation"
author: "Spencer Graves and Doug Samuelson"
date: "2020-08-04"
#output: bookdown::pdf_book
#NOTE: As of 2020-08-05 I get:
#! Package inputenc Error: Unicode character − (U+2212)
#(inputenc) not set up for use with LaTeX.
#Try other LaTeX engines instead (e.g., xelatex) if you
#are using pdflatex. See #https://bookdown.org/yihui/rmarkdown-cookbook/latex-unicode.html
# This recommends:
#output:
# pdf_document:
# latex_engine: xelatex
# 2020-09-30:
# Etienne B. Racine <notifications@github.com>
# suggested:
output:
bookdown::pdf_book:
latex_engine: xelatex
#output: bookdown::gitbook
#NOTE: As of 2019-11 "Figure \@ref(fig:plot)"
# does not work with
# output: html_document
bibliography: nuc-references.bib
vignette: >
%\VignetteIndexEntry{Forecasting nuclear proliferation}
%\VignetteIndexEntry{Time to next new nuclear-weapon state}
%\VignetteKeyword{nuclear-weapon states}
%\VignetteEngine{knitr::knitr}
%\SweaveUTF8
\usepackage[utf8](inputenc)
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# Ecdat::nuclearWeaponStates include all
# first test up to nukeDataCurrent:
nukeDataCurrent <- as.Date('2020-08-04')
library(lubridate)
library(Ecdat)
firstTstYr <- year(nuclearWeaponStates$firstTest)
(firstYear <- firstTstYr[1])
nNucStates <- nrow(nuclearWeaponStates)
nNucSt1 <- (nNucStates-1)
(Today <- today())
#currentYear <- year(Today)
currentYear <- year(nukeDataCurrent)
if((month(nukeDataCurrent)<7) ||
(difftime(nukeDataCurrent,
tail(nuclearWeaponStates$firstTest, 1),
units = 'days')<(366/2)))
currentYear <- (year(nukeDataCurrent)-1)
nYrs <- currentYear - firstYear
save_svg <- !fda::CRAN()
simStart <- (currentYear+1)
simEnd <- (currentYear+nYrs)
```
# Abstract
This article models the time between the first test of a nuclear weapon by one nation and the next over the `r nYrs` years of history since the first such test by the US available as this is being written.^[Obviously, this considers the data available as this is being written `r nukeDataCurrent`.] We use those results to forecast nuclear proliferation over the next `r nYrs` years. The maximum likelihood estimate of the time between "first tests" (`r nYrs` years divided by `r nNucSt1` "first tests") is `r round(nYrs/nNucSt1, 1)` years using the standard formula for censored estimation of a constant exponential distribution. However, a plot of the times between "first tests" of the `r nNucStates` nuclear powers as of `r nukeDataCurrent` suggests an inhomogeneous [renewal process](https://en.wikipedia.org/wiki/Renewal_theory) with the Poisson mean of the number of "first tests" each year by new nuclear-weapon states decreasing over time. This might be modeled using [`glm(..., family=poisson)`](https://www.rdocumentation.org/packages/stats/versions/3.6.1/topics/glm). Unfortunately, the standard estimate of a linear trend in log(Poisson mean) is not statistically significant. We therefore use Bayesian Model Averaging (BMA), considering two related BMA mixtures. The median of Monte Carlo simulations for the constant-linear mixture forecasts an additional 7.2 new nuclear-weapon states by `r simEnd`, using a forecasting period to match the available history,^[assuming the data used are accurate.] within the anticipated lifespan of most babies born today. Adding quadratic, cubic and quartic terms to the mixture reduces the median number of new nuclear-weapon states by `r simEnd` to 5.4. Eighty percent percent prediction limits for the quartic mixture run from 0 to 15 new nuclear-weapon states. This suggests that the risk of a nuclear war leading to the extinction of civilization is _increasing_. Nuclear proliferation will likely continue until it becomes effectively impossible for anyone to make more nuclear weapons for a very long time. This could come from a nuclear war or a massive and unprecedented strengthening of international law that provides effective judicial recourse for grievances of the poor, weak and disfranchised.
# Introduction
A plot of times between "first tests" by the world's nuclear-weapon states as of `r nukeDataCurrent` suggests that the process of nuclear proliferation has slowed; see Figure \@ref(fig:plot).
```{r plot, fig.cap = paste("Time between new nuclear-weapon states.", NucStates)}
NucStates <- paste0(
'CN = China, FR = France, GB = UK, ',
' IL = Israel,\nIN = India, KP = North Korea, ',
'PK = Pakistan, RU = Russia')
library(Ecdat)
data(nuclearWeaponStates)
ymax <- max(
nuclearWeaponStates$yearsSinceLastFirstTest,
na.rm=TRUE)
ylim0 <- c(0, ymax)
NPTdate = as.Date(c('1970-03-05', '1988-06-01'))
plotNucStates <- function(type.='n', xlim., ylim.,
line_mtext=3:2, cex.=1, mtext.=TRUE, log.='',
yNPT=NULL, ...){
##
## Write a function to create this desired plot
## that is general enough to be customized
## to make other similar but different plots
## later.
##
## Obviously, during the process of writing
## this vignette, it requires revising this
## function later as the needs become clearer.
##
## The advantage of doing it this way is that
## it makes the code easier to read, because
## it's clearer what is the same and what is
## different between similar plots.
##
# Start with an internal function
# to add the 2-letter country codes.
addCountries <- function(line_mtext=3:2, cex.=1,
mtext.=TRUE){
# Add the country codes ("ctry") to a plot
# showing the time between "first tests"
# of nuclear-weapon states
# ... to save copying code
# and hopefully make the logic clearer
xlab. <- paste(c(
'Note: The US is not on this plot,',
'because it had no predecessors.'),
collapse='\n')
if(mtext.){
mtext(xlab., 1, line_mtext[1], cex=cex.)
mtext('years from the\nprevious "first test"',
2, line_mtext[2], cex=cex.)
}
with(nuclearWeaponStates,
text(firstTest, yearsSinceLastFirstTest,
ctry, xpd=TRUE, cex=cex.))
}
# xlim and ylim?
if(missing(xlim.)){
xlim. <- range(
nuclearWeaponStates$firstTest)
}
if(missing(ylim.))ylim. <- range(
nuclearWeaponStates$yearsSinceLastFirstTest[-1])
# If very wide log scale on y,
# make the margins wider and move the label out:
if((log.=='y') && (diff(log(ylim.))>5) ){
op <- par(mar=c(5, 6, 4, 2)+0.1)
on.exit(par(op))
line_mtext[2] <- 4
}
#
plot(yearsSinceLastFirstTest~firstTest,
nuclearWeaponStates, type=type.,
xlab='', ylab='', las=1,
xlim=xlim., ylim=ylim., log=log., bty='n',
...)
abline(v=NPTdate, lty='dashed', col='grey')
if(is.null(yNPT)){
yNPT <- {
if(log.=='y') sqrt(ylim.[1]*ylim.[2]) else
mean(ylim.)
}
}
text(NPTdate-.017*diff(xlim.), yNPT,
c('NPT', 'INF'), col='grey', srt=90)
addCountries(line_mtext=line_mtext, cex.=cex.,
mtext.=mtext.)
}
plotNucStates(type.='h', ylim.=ylim0)
if(save_svg){
print(getwd())
svg('NucWeaponStates_YrsBetw1stTsts.svg')
op <- par(cex=1.7, cex.axis=1.4)
plotNucStates(type.='h', ylim.=ylim0,
mtext. = FALSE)
dev.off()
}
```
This plot also marks the effective dates of both the [Treaty on the Non-Proliferation of Nuclear Weapons (Non-Proliferation Treaty, NPT)](https://en.wikipedia.org/wiki/Treaty_on_the_Non-Proliferation_of_Nuclear_Weapons) and the [Intermediate-range Nuclear Forces (INF) Treaty](https://en.wikipedia.org/wiki/Intermediate-Range_Nuclear_Forces_Treaty), (1970-03-05 and 1988-06-01, respectively), because of the suggestion that those treaties may have slowed the rate of nuclear proliferation.
A visual analysis of this plot suggests that nuclear proliferation is continuing, and neither the NPT nor the INF treaty had a major impact on nuclear proliferation. The image is pretty bad: There were only 5 nuclear-weapon states when the NPT entered into force in 1970.^[@UN:1970. See also @Wikip:NPT.] When US President [George W. Bush](https://en.wikipedia.org/wiki/George_W._Bush) decried an ["Axis of evil"](https://en.wikipedia.org/wiki/Axis_of_evil) in his State of the Union message, 2002-01-29,^[@Bush:2002; see also @Wikip:AxisOfEvil.] there were 8. As this is written `r nukeDataCurrent`, there are 9. _Current international policy on nuclear weapons seems to assume that nuclear proliferation has stopped. It clearly has not._
@Toon:2007 noted that in 2003 another 32 countries had sufficient fissile material to make nuclear weapons if they wished. Moreover, those 32 do _NOT_ include either Turkey nor Saudi Arabia. On 2019-09-04, Turkish President Erdogan said it was unacceptable for nuclear-armed states to forbid Turkey from acquiring its own nuclear weapons.^[@Toksabay:2019; @OConnor:2019.]
Similarly, in 2006 *Forbes* reported that Saudi Arabia has "a secret underground city and dozens of underground silos for" Pakistani nuclear weapons and missiles.^[@Forbes:2006; see also @Wikip:SaudiNuc.] In 2018 the *Middle East Monitor* reported that "Israel 'is selling nuclear information' to Saudi Arabia".^[@MEM:SaudiNuc; see also @Wikip:SaudiNuc.] This is particularly disturbing, because of the substantial evidence that Saudi Arabia may have been and may still be the primary recruiter and funder of Islamic terrorism.^[@Benjamin:2016; see also @Wikiv:WinTerror.]
This analysis suggests that the number of nuclear-weapon states will likely continue to grow until some dramatic break with the past makes further nuclear proliferation either effectively impossible or sufficiently undesirable.
This vignette first reviews the data and history on this issue. We then model these data as a series of annual Poisson observations of the number of states conducting a first test of a nuclear weapon each year (1 in each of 8 years since 1945; 0 in the others).
A relatively simple model for the inhomogeneity visible in Figure \@ref(fig:plot) is [Poisson regression](https://en.wikipedia.org/wiki/Poisson_regression) assuming that log(Poisson mean) is linear in the time since the first test of a nuclear weapon by the US on `r nuclearWeaponStates$firstTest[1]`. We estimate this using [`glm(..., family=poisson)`](https://www.rdocumentation.org/packages/stats/versions/3.6.1/topics/glm). This model is plausible to the extent that this trend might represent a growing international awareness of the threat represented by nuclear weapons including a hypothesized increasing reluctance of existing nuclear-weapon states to share their technology. The current process of ratifying the new [Treaty on the Prohibition of Nuclear Weapons](https://en.wikipedia.org/wiki/Treaty_on_the_Prohibition_of_Nuclear_Weapons) supports the hypothesis of such a trend, while the lack of universal support for it and the trend visible in Figure \@ref(fig:plot) clearly indicate that nuclear proliferation is still likely to continue. We use this model to extend the `r nYrs` years of history of nuclear proliferation available as this is being written `r nukeDataCurrent` into predicting another `r nYrs` years into the future.
# How did the existing nuclear-weapon states develop this capability?
There are, of course, multiple issues in nuclear proliferation. A new nuclear-weapon state requires at least four distinct things to produce a nuclear weapon: motivation, money, knowledge, and material. The accompanying table summarizes the literature the present authors found that seems relevant to the questions at hand. Of course, definitive answers to these questions are still locked away as official secrets or have been lost to history.
However, this analysis should be sufficient to support the general conclusions of this article.
```{r howHistMat, fig.cap="Motivation, money, knowledge and material for nuclear-weapon states"}
howHist <- t(as.matrix(data.frame(
US = c('US', 'Nazi threat', 'self',
paste('own scientists + immigrants,',
'esp. fr. Germany and Italy, plus',
'collaboration w the UK & Canada'),
'Congo + self'),
"USSR" = c("USSR\n(RU)", paste(
'Hiroshima & Nagasaki bombs +',
'western invasions during WW II,',
'after WW I, and before'), 'self',
paste('own scientists + espionage',
'in the US & captured Germans'), 'self'),
"UK" = c("UK\n(GB)", 'USSR', 'self',
'Manhattan Project', 'Canada'),
"France" = c("France\n(FR)",
'USSR + Suez Crisis', rep('self', 3)),
"China" = c("China\n(CN)", paste(
'1st Taiwan Strait Crisis 1954–1955,',
'Korean Conflict, etc.'),
'self', 'USSR', 'self'),
"India" = c("India\n(IN)",
paste('loss of territory:',
'China-Himalayan border-1962'),
'self',
'students in UK, US',
'Canadian nuc reactor'),
"Israel" = c("Israel\n(IL)",
'hostile neighbors', 'self',
'self + France', 'France + ???'),
"Pakistan" =c("Pakistan\n(PK)",
'Loss of E. Pakistan in 1971',
'Saudis + self', 'US, maybe China?', 'self?'),
"N.Korea" = c("N.Korea\n(KP)",
'threats fr. US', "self?",
'US via Pakistan?', 'self?'),
stringsAsFactors=FALSE)))
colnames(howHist) <- c('Nation', 'Motivation',
'Money', "Knowledge", "Material")
library(knitr)
library(kableExtra)
emTot <- 41
em. <- c(nation=2.2, motive=10, money=2,
knowledge=10, material=3.4)
em <- paste0(emTot*em./sum(em.), 'em')
kable(howHist, row.names=FALSE) %>%
column_spec(1, width=em[1]) %>%
column_spec(2, width=em[2]) %>%
column_spec(3, width=em[3]) %>%
column_spec(4, width=em[4]) %>%
column_spec(5, width=em[5])
```
## Motivation
Virtually any country that feels threatened would like to have some counterweight against aggression by a potential enemy.
- The US funded the Manhattan project believing that Nazi Germany likely had a similar project.
- Soviet leaders might have felt a need to defend themselves from nuclear coercion after having been invaded by Nazi Germany only a few years earlier, and having defeated [foreign invasions from the West and the East after world War I trying to put the Tsar back in power](https://en.wikipedia.org/wiki/Allied_intervention_in_the_Russian_Civil_War).^[@Fogelsong:1995. That doesn't count [numerous other invasions that are a sordid part of Russian history](https://en.wikipedia.org/wiki/French_invasion_of_Russia), which educated Russians throughout history would likely remember, even if their invaders may not.]
- The United Kingdom and France felt nuclear threats from the Soviet Union.^[The UK and France would have had many reasons to fear the intentions of the USSR during the early period of the [Cold War](https://en.wikipedia.org/wiki/Cold_War): The first test of a nuclear weapon by the USSR came just over three months after the end of the 1948-49 [Berlin Blockade](https://en.wikipedia.org/wiki/Berlin_Blockade). Other aspects of Soviet repression in countries they occupied in Eastern Europe contributed to the failed [Hungarian Revolution of 1956](https://en.wikipedia.org/wiki/Hungarian_Revolution_of_1956).]
- France's concern about the Soviets increased [after the US refused to support them](https://en.wikipedia.org/wiki/France_and_weapons_of_mass_destruction#cite_note-16) during the [1956 Suez Crisis](https://en.wikipedia.org/wiki/Suez_Crisis): If the US would not support a British-French-Israeli invasion of Egypt, the US might not defend France against a possible Soviet invasion.^[@Fromkin:2006. See also @Wikip:FranceNuc]
- China reportedly decided to initiate its nuclear weapons program during the [First Taiwan Strait Crisis](https://en.wikipedia.org/wiki/China_and_weapons_of_mass_destruction#Nuclear_weapons) of 1954-55,^[@Wikip:ChinaNuc; see also @Wikip:TaiwanStrait, @Halperin:1966, and @Wikip:Ellsberg.] following nuclear threats from the US regarding Korea.^[@Pierson:2017. See also @Wikip:Ellsberg.]
- India lost territory to China in the [1962 Sino-Indian War](https://en.wikipedia.org/wiki/Sino-Indian_War), which reportedly convinced India to abandon a policy of avoiding nuclear weapons.^[@Riedel:2012. See also @Wikip:IndiaNuc. India and China have continued to have conflicts. See, for example, the Wikipedia articles on [China-India relations](https://en.wikipedia.org/wiki/China%E2%80%93India_relations) and the [2017 China-India border standoff](https://en.wikipedia.org/wiki/2017_China%E2%80%93India_border_standoff).]
- Pakistan's nuclear weapons program began in 1972 in response to the loss of East Pakistan (now Bangladesh) in the 1971 [Bangladesh Liberation War](https://en.wikipedia.org/wiki/Bangladesh_Liberation_War).^[@Wikip:PakNuc. [India-Pakistan relations](https://en.wikipedia.org/wiki/India%E2%80%93Pakistan_relations) have been marked by frequent conflict since the two nations were born with the dissolution of the British Raj in 1947. This history might help people understand the need that Pakistani leaders may have felt and still feel for nuclear parity with India, beyond the loss of half their population and 15 percent of their land area in the 1971 Bangladesh Liberation War.] On November 29, 2016, Moeed Yusuf claimed that the threat of a nuclear war between India and Pakistan was the most serious foreign policy issue facing then-President-elect Trump.^[@Yusuf:2016] That may have been an overstatement, but the possibilities of a nuclear war between India and Pakistan should not be underestimated. [There have been lethal conflicts between India and Pakistan at least as recent as 2019.](https://en.wikipedia.org/wiki/Indo-Pakistani_wars_and_conflicts#Past_skirmishes_and_standoffs) If that conflict goes nuclear, it could produce a “nuclear autumn” during which a quarter of humanity not directly impacted by the nuclear war would starve to death, according to simulations by leading climatologists.^[@Toon:2007]
- Israel has faced potentially hostile neighbors since its declaration of independence in 1948.^[@Wikip:ArabIsrael. Threats perceived by Israel continue, including the [Gaza border protests](https://en.wikipedia.org/wiki/2018_Gaza_border_protests) that have continued at least into 2020. One might therefore reasonably understand why Israel might feel a need for nuclear weapons and why others might believe that the 1979-09-22 [Vela incident](https://en.wikipedia.org/wiki/Vela_Incident) was an Israeli nuclear test.]
- North Korea first tested a nuclear weapon on 2006-10-09,^[@CRS:2016. The US Congressional Research Service in 2016 reported, “The Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) [`PrepCom`](https://en.wikipedia.org/wiki/Preparatory_Commission_for_the_Comprehensive_Nuclear-Test-Ban_Treaty_Organization)'s international monitoring system detected data indicating that North Korea had conducted a nuclear test on January 6, 2016. ... On October 9, 2006, North Korea declared that it had conducted an underground nuclear test.” For the present purposes, we use the October date declared by North Korea, not the January date reported by CTBTO. See also @Wikip:NKoreaNuc06.] less than five years after having been named as part of an ["Axis of evil"](https://en.wikipedia.org/wiki/Axis_of_evil) by US President George W. Bush on 2002-01-29.^[@Bush:2002; see also @Wikip:AxisOfEvil.] Chomsky claimed that the relations between the US and North Korea have followed "a kind of tit-for-tat policy. You make a hostile gesture, and we'll respond with a crazy gesture of our own. You make an accommodating gesture, and we'll reciprocate in some way." He gave several examples including a 1994 agreement that halted North Korean nuclear-weapons development. "When George W. Bush came into office, North Korea had maybe one [untested] nuclear weapon and verifiably wasn't producing any more."^[@Chomsky:2017, pp. 131-134. Chomsky includes in this game of tit-for-tat the total destruction of North Korean infrastructure during the Korean War in the early 1950s, including huge dams that controlled the nation's water supply, destroying their crops, and raising the specter of mass starvation. @Kolko:1968 noted that German General [Syss-Inquart](https://en.wikipedia.org/wiki/Arthur_Seyss-Inquart) ordered similar destruction of dikes in Holland in 1945, which condemned many Dutch civilians to death by starvation. For that crime Syss-Inquart became one of only 24 of the people convicted at the Nuremberg war crimes trial to have been sentenced to death. Chomsky noted that this is "not in our memory bank, but it's in theirs."]
All this suggests that it will be difficult to reduce the threat of nuclear proliferation and nuclear war without somehow changing the nature of international relations so weaker countries have less to fear from the demands of stronger countries.
## Money
To help us understand the differences in sizes of the different nuclear-weapon states, Figure \@ref(fig:nucStates) plots the populations and Gross Domestic Products (GDP) of the current nuclear-weapon states.^[Data for different years, 2017-2020, depending on what was available from Wikipedia on 2020-02-05.] The following subsections provide analysis with references behind the summary in Figure \@ref(fig:nucStates).
```{r nucStates, fig.cap = "Gross Domestic Product and Population of Nuclear-Weapon States. (Country codes as with Figure 1 with GB, FR overplotted.)"}
#op <- par(mar=c(5,5,4,2)+.1)
plot(GDP_B/1000 ~ popM, nuclearWeaponStates, type='n',
log='xy', las=1, xlab='', ylab='')
title(ylab='nominal GDP (USD trillions)')
title(xlab='population (millions)', line = 1.7)
nNucStates <- nrow(nuclearWeaponStates)
i0 <- c(1:2, 4:(nNucStates-1))
with(nuclearWeaponStates, arrows(
popM[i0], GDP_B[i0]/1000,
popM[i0+1], GDP_B[i0+1]/1000, col='grey',
lty='dotted'))
cols <- c(US='blue', RU='red', GB='red',
FR='blue', CN='red', IN='orange',
IL='blue', PK='green', KP='red')
with(nuclearWeaponStates,
text(popM, GDP_B/1000, ctry,
col=cols) )
leg <- with(nuclearWeaponStates,
paste(ctry, '=', nation))
#with(nuclearWeaponStates,
# legend('bottomright', legend=leg,
# bty='n', text.col=cols))
#par(op)
if(save_svg){
print(getwd())
svg('NucWeaponStates_GDP_pop.svg')
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,4, 2, 2)+.1)
cex2 <- 1.7
plot(GDP_B/1000 ~ popM, nuclearWeaponStates, type='n',
log='xy', las=1, xlab='', ylab='')
nNucStates <- nrow(nuclearWeaponStates)
i0 <- c(1:2, 4:(nNucStates-1))
with(nuclearWeaponStates, arrows(
popM[i0], GDP_B[i0]/1000,
popM[i0+1], GDP_B[i0+1]/1000, col='grey',
lty='dotted'))
cols <- c(US='blue', RU='red', GB='red',
FR='blue', CN='red', IN='orange',
IL='blue', PK='green', KP='red')
with(nuclearWeaponStates,
text(popM, GDP_B/1000, ctry,
col=cols) )
dev.off()
}
```
It's no accident that most of the world's nuclear-weapon states are large countries with substantial populations and economies. That's not true of Israel with only roughly 9 million people nor North Korea with roughly 26 million people in 2018. France and the UK have only about 67 and 68 million people, but they are also among the world leaders in the size of their economies.
Pakistan is a relatively poor country. It reportedly received financial assistance from Saudi Arabia for its nuclear program.^[@Riedel:2008.]
Another reason for a possible decline in the rate of nuclear proliferation apparent in Figure \@ref(fig:plot) is the fact that among nuclear-weapon states, those with higher GDPs tended to acquire this capability earlier, as is evident in Figure \@ref(fig:nucStates).
## Knowledge
In 1976, [John Aristotle Phillips](https://en.wikipedia.org/wiki/John_Aristotle_Phillips), an "underachieving" undergraduate at Princeton University, "designed a nuclear weapon using publicly available books and papers."^[@Spokane:1976. See also @Wikip:studentBomb.] Nuclear weapons experts disagreed on whether the design would have worked. Whether Phillips' design would have worked or not, it should be clear that the continuing progress in human understanding of [nuclear physics](https://en.wikipedia.org/wiki/Nuclear_physics) inevitably makes it easier for people interested in making such weapons to acquire the knowledge of how to do so.
Before that, the nuclear age arguably began with the 1896 discovery of radioactivity by the French scientist Henri Becquerel. It was further developed by Pierre and Marie Curie in France, Ernest Rutherford in England, and others, especially in France, England and Germany.^[@Wikip:NucPhys] In 1933 after Adolf Hitler came to power in Germany, Leo Szilard moved from Germany to England. The next year he patented the idea of a nuclear fission reactor. After World War II began, the famous [Manhattan Project](https://en.wikipedia.org/wiki/Manhattan_Project) became a joint British-American project, which produced the very first test of a nuclear weapon.^[@Wikip:HistNucWeapons]
After Soviet premier Joseph Stalin learned of the atomic bombings of Hiroshima and Nagasaki, the USSR (now Russia) increased the funding for their nuclear-weapons program. That program was helped by intelligence gathering about the German nuclear weapons project and the American Manhattan Project.^[@osti:espionage. See also @wikip:Soviet.]
The UK's nuclear-weapons program was built in part on their wartime participation in the Manhattan Project, as noted above. France was among the leaders in nuclear research until World War II. They still had people with the expertise needed after the Suez Crisis convinced them they needed to build nuclear bombs, as noted above.^[See also @Wikip:HistNucWeapons]
China got some help from the Soviet Union during the initial phases of their nuclear program.^[@Wikip:ChinaNuc.]
The first country to get nuclear weapons after the Non-Proliferation Treaty was India. Their Atomic Energy Commission was founded in 1948, chaired by Homi J. Bhabha. He had published important research in nuclear physics while a graduate student in England in the 1930s, working with some of the leading nuclear physicists of that day.^[@Wikip:Bhabha; see also @Wikip:NucTimeline.]
Meanwhile, Israel's nuclear weapons program initially included sending students abroad to study under leading physicists like Enrico Fermi at the University of Chicago. It also included extensive collaboration with the French nuclear-weapons program.^[@Wikip:IsraelNuc.]
Pakistan got [secret help from the US in the 1980s in violation of US law to secure Pakistani cooperation with US support for anti-Soviet resistance in Afghanistan](https://www.wilsoncenter.org/publication/new-documents-spotlight-reagan-era-tensions-over-pakistani-nuclear-program).[^Pak2] [Robert Gallucci](https://en.wikipedia.org/wiki/Robert_Gallucci) said that the nuclear weapons programs of North Korea, Iran and Libya would not have gotten off the ground without help the US secretly gave to the Pakistani nuclear weapons program in violation of US law; Gallucci is a leading scholar and expert on non-proliferation, who held senior positions in the Reagan, George H. W. Bush and Clinton administrations.[^KP] [Vikram Sood](https://en.wikipedia.org/wiki/Vikram_Sood), a former head of India's foreign intelligence agency, echoed Galluci's claims, adding that Pakistan ''may'' have given nuclear-weapons technology to al Qaeda "just weeks prior to September 11, 2001."^[Sood:2008] It may not be wise to accept Sood's claim at face value, given the long-standing hostility between India and Pakistan. On the other hand, we should not dismiss the possibility that a terrorist organization might acquire nuclear weapons.
Western sources have claimed that China also helped Pakistan's nuclear-weapons program, but China has denied those claims.^[@Wikip:PakNuc.]
And now the US is helping Saudi Arabia obtain nuclear power, in spite of (a) the evidence that [the Saudi government including members of the Saudi royal family were involved at least as early as 1999 in preparations for the suicide mass murders of September 11, 2001](https://en.wikipedia.org/wiki/The_28_pages),^[@Graham:2003. See also @Wikip:28pages.] and (b) their [on-going support for Al Qaeda in Yemen, reported as recently as 2019](https://en.wikipedia.org/wiki/Saudi_Arabian-led_intervention_in_Yemen).^[See, for example, @Bazzi:2019 and @Wikip:SaudiYemen, more generally.]
## Material
Reportedly the most difficult part of making nuclear weapons today is obtaining sufficient fissile material. @Toon:2007 said, "Thirteen countries operate plutonium and/or uranium enrichment facilities, including Iran", but Iran did not have sufficient fissile material in 2003 to make a nuclear weapon. Another 20 were estimated to have had sufficient stockpiles of fissile material acquired elsewhere to make nuclear weapons. They concluded that 32 (being 13 minus 1 plus 20) additional countries have sufficient fissile material to make nuclear weapons if they want.^[pp. 1975, 1977. The 32 countries they identified included 12 of the 13 that "operate plutonium and/or uranium enrichment facilities", excepting Iran as noted. The other 20 countries acquired stockpiles elsewhere. In addition to the 32 with sufficient fissile material to make a nuclear weapon, Egypt, Iraq and the former Yugoslavia were listed as having abandoned a nuclear-weapons program.]
@Toon:2007 also said, "In 1992 the International Atomic Energy Agency safeguarded less than 1% of the world’s HEU [Highly Enriched Uranium] and only about 35% of the world inventory of Pu [Plutonium] ... . Today [in 2007] a similarly small fraction is safeguarded."
HEU is obtained by separating $^{235}$U, which is only 0.72 percent of naturally occurring uranium.^[@GlobalSec:Uriso.] Weapons-grade uranium has at least 85 percent $^{235}$U.^[@Wikip:EnrichedU, section on "Highly enriched uranium (HEU)".] Thus, at least 0.85/0.0072 = 118 kg of naturally occurring uranium are required to obtain 1 kg that is weapons-grade, and @Toon:2007 estimated that 25 kg of HEU would be used on average for each $^{235}$U-based nuclear weapon. Plutonium, by contrast, is a byproduct of energy production in standard $^{238}$U nuclear reactors.
Much of the uranium for the very first test of a nuclear weapon by the US came from the Congo,^[@Wikip:Manhattan.] but domestic sources provided most of the uranium for later US nuclear-weapons production.^[@Wikip:Ures.] The Soviet Union (USSR, now Russia) also seems to have had adequate domestic sources for its nuclear-weapons program, especially including Kazakhstan, which was part of the USSR until 1990; Kazakhstan has historically been the third largest source of uranium worldwide after Canada and the US.^[@Wikip:Ures.] The UK presumably got most of its uranium from Canada.
The French nuclear-weapons program seems to have been built primarily on plutonium.^[@Wikip:FranceNuc. See also Table 2 in @Toon:2007, which claims that in 2003, France had enough fissile material for roughly 24,000 plutonium bombs and 1,350 $^{235}$U bombs.] This required them to first build standard $^{238}$U nuclear reactors to make the plutonium. Then they didn't need nearly as much uranium to sustain their program.
China has reportedly had sufficient domestic reserves of uranium to support its own needs,^[@Wikip:Ures.] even exporting some to the USSR in the 1950s in exchange for other assistance with their nuclear defense program.^[@Wikip:ChinaNuc.]
India's nuclear weapons program seems to have been entirely (or almost entirely) based on plutonium.^[@Wikip:IndiaNuc; see also @Toon:2007 and @Wikip:Ures.]
Israel seems not to have had sufficient uranium deposits to meet its own needs. Instead, they purchased some from France until France ended their nuclear-weapons collaboration with Israel in the 1960s. To minimize the amount of uranium needed, nearly all Israeli nuclear weapons seem to be plutonium bombs.^[@Toon:2007.]
It's not clear where Pakistan got most of its uranium: Its reserves in 2015 were estimated at zero, and its historical production to that point was relatively low.^[@Wikip:Ures.] By comparison with the first seven nuclear-weapon states, it's not clear where Pakistan might have gotten enough uranium to produce 83 plutonium bombs and 44 uranium bombs, as estimated by @Toon:2007^[Table 2, p. 1976.] As previously noted, the US helped the Pakistani nuclear-weapons program in the 1980s and accused China of providing similar assistance, a charge that China has repeatedly and vigorously denied. China has provided civilian nuclear reactors, which could help produce plutonium but not $^{235}$U.^[@Wikip:PakNuc.]
According to the Federation of American Scientists, "North Korea maintains uranium mines with an estimated four million tons of exploitable high-quality uranium ore ... that ... contains approximately 0.8% extractable uranium."^[@Aftergood:2006; see also @Wikip:NKoreaNuc.] If that's accurate, processing all that would produce 4,000,000 times 0.008 = 32,000 tons of pure natural uranium, which should be enough to produce the weapons they have today.
## Conclusions regarding motivation, money, knowledge and material
1. There seems to be no shortage of motivations for other countries to acquire nuclear weapons. The leaders of the Soviet Union had personal memories of being invaded not only by Germany during World War II but also by the US and others after World War I. The UK had reason to fear the Soviets in their occupation of Eastern Europe. The French decided after Suez they couldn't trust the US to defend them. China had been forced to yield to nuclear threats before starting their nuclear program, as did India, Pakistan and North Korea. Israel has fought multiple wars since their independence in 1948.
2. The knowledge and material required to make such weapons in a relatively short order are also fairly widely available, even without the documented willingness of current nuclear powers to secretly help other countries acquire such weapons in some cases.[^quickNukes]
3. Unless there is some fundamental change in the structure of international relations, it seems unwise to assume that there will not be more nuclear-weapon states in the future, with the time to the next "first test" of a nuclear weapon following a probability distribution consistent with the previous times between "first tests" of nuclear weapons by the current nuclear-weapon states.
# Distribution of the time between Poisson "first tests"
Possibly the simplest model for something like the time between "first tests" in an application like this is to assume they come from one [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution) with 8 observed times between the 9 current nuclear-weapon states plus one [censored observation](https://en.wikipedia.org/wiki/Censoring_(statistics)) of the time between the most recent one and a presumed next one. This simple theory tells us that the maximum likelihood estimate of the mean time between such "first tests" is the total time from the US "Trinity" test to the present, `r round(yrsSinceT <- as.numeric(difftime(nukeDataCurrent, nuclearWeaponStates$firstTest[1], units='days')/365.25), 1)` years, divided by the number of new nuclear-weapon states, `r (nNewNucSt <- nNucStates-1)`, not counting the first, which had no predecessors. Conclusion: Mean time between "first tests" = `r round(yrsSinceT / nNewNucSt, 1)` years.[^WikivTime]
However, Figure \@ref(fig:plot) suggests that the time between "first tests" of succeeding nuclear-weapon states is increasing. The decreasing hazard suggested by this figure requires mathematics that are not as easy as the censored data estimation as just described.
To understand the current data better, we redo Figure \@ref(fig:plot) with a log scale on the y axis in Figure \@ref(fig:ploty).
```{r ploty, fig.cap = paste("Semilog plot of time between new nuclear-weapon states.", NucStates)}
plotNucStates(log.='y')
# optionally write to a file
if(save_svg){
print(getwd())
svg('NucWeaponStates_logYrsBetw1stTsts.svg')
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 2)+.1)
cex2 <- 1.7
plotNucStates(mtext. = FALSE, log.='y',
cex.=cex2)
dev.off()
}
```
Figures \@ref(fig:plot) and \@ref(fig:ploty) seem consistent with the following:
- If the mean time between "first tests" is increasing over time, as suggested by Figures \@ref(fig:plot) and \@ref(fig:ploty), then the distribution cannot be exponential, because that requires a constant [hazard rate](https://en.wikipedia.org/wiki/Survival_analysis#Hazard_function_and_cumulative_hazard_function). [For the exponential distribution,
$$h(t) = (-d/dt \log S(t)) = \lambda$$,
writing the exponential survival function as $S(t) = \exp(-\lambda t)$.]
- Even though nuclear proliferation has been slowing since 1950, it could _accelerate_ in the future if more states began to perceive greater threats from other nations.
- Fortunately we can simplify this modeling problem by using the famous duality between exponential time between events and a Poisson distribution for numbers of events in specific intervals of time. By modeling Poisson counts of "first tests" each year, we can use techniques for Poisson regression for models suggested by Figure \@ref(fig:ploty). The simplest such model might consider log(Poisson mean numbers of "first tests" each year) to be linear in the time since the first test of a nuclear weapon (code-named ["Trinity"](https://en.wikipedia.org/wiki/Trinity_(nuclear_test))).^[@Rhodes:1986. See also @Wikip:Trinity.] However, the image in Figure \@ref(fig:ploty) suggests the line may not be straight. Easily tested alternatives to linearity could be second, third and fourth powers of the `timeSinceTrinity`.^[One might also consider a model with the log(Poisson mean) behaving like a ["Wiener process" (also called a "Brownian motion")](https://en.wikipedia.org/wiki/Wiener_process). This stochastic formulation would mean that the variance of the increments in log(hazard) between "first tests" is proportional to the elapsed time. See @Wolfram:Wiener and @Wikip:Wiener. The [`bssm` package](https://www.rdocumentation.org/packages/bssm/versions/0.1.7) provides a reasonable framework for modeling this. Its [`ng_bsm` function](https://www.rdocumentation.org/packages/bssm/versions/0.1.7/topics/ng_bsm) supports modeling a normal random walk in log(Poisson mean) of the number of "first tests" each year. In this article, we model the trend as deterministic and leave consideration of a [Gaussian random walk](https://en.wikipedia.org/wiki/Random_walk) and similar stochastic formulations for future work.]
We use Poisson regression to model this as a series of the number of events each year.^[We could potentially use one observation each month, week or day. Such a change might give us slightly better answers while possibly increasing the compute time more than it's worth.]
# Parameter estimation
For modeling and parameter estimation, we use [`glm(firstTests ~ timeSinceTrinity, poisson)`](https://www.rdocumentation.org/packages/stats/versions/3.6.1/topics/glm) with:
- `firstTests` = the number of "first tests" of a nuclear-weapon by a new nuclear-weapon state each year, and
- `timeSinceTrinity` = number of years since 1945-07-16, when the first nuclear weapon was tested, code-named ["Trinity"](https://en.wikipedia.org/wiki/Trinity_(nuclear_test)).
We use the [`lubridate` package](https://lubridate.tidyverse.org/) for dates. The first thing we want is the most recent date we can confidently use for the validity of the dataset we use, `Ecdat::nuclearWeaponStates`. This is in a `Date` variable `nukeDataCurrent` defined as `r nukeDataCurrent` in a `setup` section before the visible start of the narrative of this vignette. From this we get the year:
```{r year}
(currentYear <- year(nukeDataCurrent))
```
We include an observation for the current year only if it's more than 6 months since January 1 and since the last "first test".
```{r currentYear}
if((month(nukeDataCurrent)<7) ||
(difftime(nukeDataCurrent,
tail(nuclearWeaponStates$firstTest, 1),
units = 'days')<(366/2)))
currentYear <- (year(nuclearWeaponStates)-1)
```
Start after the year of the first test of a nuclear weapon.
```{r firstYear}
firstTstYr <- year(nuclearWeaponStates$firstTest)
(firstYear <- firstTstYr[1])
```
We use this to create a vector of the number of `firstTest`s by year and put this in a `tibble` with `Year`.
```{r firstTests}
(nYrs <- currentYear - firstYear)
firstTests <- ts(rep(0, nYrs), firstYear+1)
firstTstYrSinceFirst <- firstTstYr - firstYear
firstTests[firstTstYrSinceFirst] <- 1
library(tibble)
(FirstTsts <- tibble(Year=time(firstTests),
nFirstTests=firstTests))
```
We add `ctry` to this `tibble` for future reference.
```{r firstCtry}
Ctry <- rep('', nYrs)
Ctry[firstTstYrSinceFirst] <-
nuclearWeaponStates$ctry[-1]
FirstTests <- cbind(FirstTsts, ctry=Ctry)
```
We add `timeSinceTrinity`, which we will use in modeling.
```{r timeSince}
FirstTests$timeSinceTrinity <- 1:nYrs
```
We then fit a model with log(Poisson mean number of first tests each year) linear in `timeSinceTrinity`.
```{r fit1}
summary(fitProlif1 <- glm(
firstTests ~ timeSinceTrinity,
poisson, FirstTests))
```
This says that the time trend visible in Figures \@ref(fig:plot) and \@ref(fig:ploty) is not statistically significant.
[George Box](https://en.wikipedia.org/wiki/George_E._P._Box) famously said that, ["All models are wrong, but some are useful."](https://en.wikipedia.org/wiki/All_models_are_wrong)^[@BoxDraper:1987; @Wikip:ModelsWrong.]
@Burnham:1998 and others claim that better predictions can generally be obtained using Bayesian Model Averaging.^[See also @Raftery:1995 and @Claeskens:2008.] In this case, we have two models: `log(Poisson mean)` being constant or linear in `timeSinceTrinity`. The [`bic.glm`](http://finzi.psych.upenn.edu/R/library/BMA/html/bic.glm.html) function in the [`BMA`](https://CRAN.R-project.org/package=BMA) package can estimate these two models and compute posterior probabilities.
```{r BMA}
library(BMA)
fitProlif <- bic.glm(
FirstTests['timeSinceTrinity'],
FirstTests$nFirstTests,
"poisson")
summary(fitProlif)
```
It is standard in the BMA literature to assume a priori an approximate uniform distribution over all models considered with a penalty for estimating each additional parameter to correct for the tendency of the models to overfit the data. With these standard assumptions, this comparison of these two models estimates a `r round(100*fitProlif$postprob[2])` percent posterior probability for the model linear in `timeSinceTrinity`, leaving `r round(100*fitProlif$postprob[1])` percent probability for the model with a constant Poisson mean. Figure \@ref(fig:plotyf) adds these lines to Figure \@ref(fig:ploty).^[For Figure \@ref(fig:plotyf), we use the standard duality between the Poisson and exponential distributions. Of course, when the hazard rate is not constant, the distribution of the time to the next "first test" is not exponential. Modeling log(Poisson mean) gives us more flexibility than trying to use any of the standard generalizations of the exponential distribution.]
```{r plotyf, fig.cap = paste("BMA fit to time between new nuclear-weapon states.", NucStates)}
plotNucStates(log.='y', yNPT=4)
predProlif <- with(fitProlif,
outer(rep(1, nYrs+1), mle[, 1]) +
outer(0:nYrs, mle[, 2]))
lgnd <- paste0(c('constant', 'linear'),
' (', 100*round(fitProlif$postprob, 2), '%)')
firstTest_nYrs <- as.Date(paste0(
trunc(nuclearWeaponStates$firstTestYr[1])+0:nYrs,
'-07-01') )
matlines(firstTest_nYrs, exp(-predProlif),
lty=c('dashed', 'dotted'),
col=c('red', 'blue'))
legend('topleft', lty=c('dashed', 'dotted'),
col=c('red', 'blue'), lgnd,
bty='n')
```
The lines in this figure seem higher than the mean of the points and a linear trend through the points. This bias might be explained by the difference between ordinary least squares and `glm` used in this case.
It's well known that extrapolation is problematic. Bayesian Model Averaging offers on average more plausible predictions than using a single model. Before proceeding, let's consider a similar BMA fit with quadratic, cubic, and quartic terms.^[Predictions that do not rely on linear or polynomial extrapolation might be obtained from a Gaussian random walk model of the log(Poisson mean). The `bssm` package claims to be able to support that kind of model. However, we have so far been unable to get sensible answers from that software and have therefore omitted it from the present discussion. We doubt if it would change the overall conclusions, though some of the details would change.]
Let's first try a model with linear, quadratic, cubic, quartic and quintic terms:
```{r fit2345}
FirstTests$time2 <- (1:nYrs)^2
FirstTests$time3 <- (1:nYrs)^3
FirstTests$time4 <- (1:nYrs)^4
FirstTests$time5 <- (1:nYrs)^5
fitProlif5 <- try(bic.glm(
FirstTests[4:8], FirstTests$nFirstTests, "poisson"))
```
Evidently, `bic.glm` cannot estimate a quintic model.
```{r fit23}
fitProlif4 <- bic.glm(
FirstTests[4:7], FirstTests$nFirstTests, "poisson")
fitProlif4$postprob
fitProlif4$mle
```
When quadratic, cubic, and quartic terms are considered, the `BMA:::bic.glm` algorithm keeps only the highest order term for each model: All the models retained have an intercept and a power of 1, 2, 3 or 4 of `timeSinceTrinity`. Moreover, their regression coefficients are all negative. This means that for each model in the Poisson mixture, the minimum of the mean time to the next "first test" occurs when `timeSinceTrinity` is zero. When a fifth order term is included, one of the models the algorithm tries to fit is computationally singular. Both these results make some sense, as there are only 8 years with one "first test"; all the others have zero "first tests", and no year had more than one.
For `nYrs` = 74 and 75, only the highest order terms in each model (quadratic, cubic, quartic) were retained. Let's confirm that.
```{r highOrder}
Order <- rep(0, 5)
names(Order) <- c('const.', 'linear',
'quadratic', 'cubic', 'quartic')
if(nrow(fitProlif4$mle)!=5){
print(fitProlif4$mle)
stop('fitProlif4 does not retain exactly ',
'5 models, unlike with nYrs=74, 75...???')
}
for(pwr in 2:5){
ip <- which(fitProlif4$mle[, pwr] !=0)
if(length(ip) != 1){
print(fitProlif4$mle)
stop('power ', pwr, ' does not appear ',
'in exactly one model, ',
'unlike with nYrs=74, 75...???')
}
Order[pwr] <- (ip-1)
}
Order
```
We add the extra lines of `fitProlif4` to Figure \@ref(fig:plotyf) to get Figure \@ref(fig:ploty4).
```{r ploty4, fig.cap = paste("BMA quartic fit to time between new nuclear-weapon states.", NucStates)}
plotNucStates4 <- function(yNPT.=4, ...){
plotNucStates(log.='y', yNPT=yNPT., ...)
predProlif4 <- matrix(NA, nYrs+1, 5)
# timeSinceTrinity=0
predProlif4[1, ] <- fitProlif4$mle[,1]
# constant model
predProlif4[, 1] <- fitProlif4$mle[1,1]
colnames(predProlif4) <- names(Order)[
1+Order]
# earlier code; obsolete
#for(pwr in 1:4){
# predProlif4[, pwr+1] <- with(fitProlif4,
# mle[pwr+1, 1] + ((0:nYrs)^pwr)*
# fitProlif4$mle[pwr+1, pwr+1])
#}
predProlif4[-1, 2:5] <- (
outer(rep(1, nYrs), fitProlif4$mle[-1, 1]) +
as.matrix(FirstTests[4:7]) %*%
t(fitProlif4$mle[-1, -1]) )
lgnd4 <- paste0(c('constant', 'linear',
'quadratic', 'cubic', 'quartic'), ' (',
100*round(fitProlif4$postprob[1+Order],
4), '%)')
matlines(firstTest_nYrs,
exp(-predProlif4[, Order+1]),
lty=1:5, col=1:5)
legend('topleft', lty=1:5, col=1:5, lgnd4,
bty='n', cex=0.95)
}
plotNucStates4()
if(save_svg){
print(getwd())
svg('NucWeaponStates_BMAyrsBetw1stTsts.svg',
height=3.5)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1,
mfrow=c(1, 2))
# const+linear
xlim4 <- c(min(nuclearWeaponStates$firstTest),
max(nuclearWeaponStates$firstTest)+365)
plotNucStates(log.='y', yNPT=4, mtext. = FALSE,
xlim.=xlim4)
matlines(firstTest_nYrs, exp(-predProlif),
lty=c('dashed', 'dotted'),
col=c('red', 'blue'))
legend('topleft', lty=c('dashed', 'dotted'),
col=c('red', 'blue'), lgnd,
bty='n')
# const+...+quartic
plotNucStates4(mtext.=FALSE, xlim.=xlim4)
# plotNucStates(log.='y', yNPT=4, mtext. = FALSE)
# matlines(firstTest_nYrs, exp(-predProlif4),
# lty=1:5, col=1:5)
# legend('topleft', lty=1:5, col=1:5, lgnd4,
# bty='n', cex=0.85)
dev.off()
}
```
Comparing predictions between `fitProlif` and `fitProlif4` might help us understand better the limits of what we can learn from the available data. A visual analysis of Figure \@ref(fig:ploty4) makes one wonder if the quartic, cubic and quadratic fits are really almost as good as the linear, as suggested by minor differences in the posterior probabilities estimated by the `bic.glm` algorithm. However, the forecasts of nuclear proliferation will be dominated by the constant component of the BMA mixture; its posterior probability is `r 100*round(fitProlif$postprob[1], 2)` percent for the constant-linear mixture and `r 100*round(fitProlif4$postprob[1], 4)` percent for the quartic mixture. That means that the median line and all the lower quantiles of all simulated futures based on these models would be dominated by that constant term.
Moreover, the quadratic, cubic and quartic lines in the right (quartic mixture) panel of Figure \@ref(fig:ploty4) do not look nearly as plausible, at least to the present author, as the constant and linear lines.^[Recall that the estimation methodology here is Poisson regression, not ordinary least squares.] That, in turn, suggests that the constant linear mixture may be more plausible than the quartic mixture
We next use `fitProlif` and `fitProlif4` to compute central 60 and 80 percent confidence limits plus 80 percent prediction, and (0.8, 0.8) tolerance limits for future nuclear proliferation, as discussed in the next three sections of this vignette.[^Confidence_Interval]
# Confidence limits
We start by computing `nSims` simulated Poisson mean numbers of "first tests" by new nuclear-weapon states for each of the `nYrs` years used in `fitProlif` and `fitProlif4` and another `nYrs` years beyond. These simulations will later be used to compute confidence limits for the model estimates of the Poisson mean and prediction and tolerance limits for the actual number of nuclear-weapon states.
<!--I wanted to put "[^Confidence_Interval]" here. However, that duplicated a relatively large footnote rather than creating a second reference to the same large footnote; unacceptable.-->
```{r meanSims}
nSims <- 5000
timeSncT <- 1:(2*nYrs)
pastfut <- tibble(Year=firstYear+timeSncT,
timeSinceTrinity=timeSncT,
time2=timeSncT^2, time3=timeSncT^3,
time4=timeSncT^4)
library(Ecfun)
simMeans <- simulate(fitProlif, nSims, seed=3,
newdata=pastfut[2], type='response')
dim(simMeans)
# earlier simulations showed a curve
# for the mean of simMeans4 that
# that was quite extreme.
# Check it by set.seed(1) here
# and 2, 3, ... later
simMeans4 <- simulate(fitProlif4, nSims, seed=1,
newdata=pastfut[2:5], type='response')
dim(simMeans4)
```
We invert these simulated Poisson means to get simulated exponential times, then summarize them in a format compatible with `yearsSinceLastFirstTest` in `nuclearWeaponStates`.
```{r expSum}
sumSims <- function(x, Year=pastfut$Year){
##
## return data.frame of Year with
## mean and (.1, .5, .9) quantiles of x
##
Yr <- as.Date(paste0(Year, '-07-01'))
xMean <- apply(x, 1, mean)
xCI <- apply(x, 1, quantile,
probs=c(.1, .2, .5, .8, .9))
# fix names
rownames(xCI) <- c(
'L10', 'L20', 'median', 'U20', 'U10')
xSum <- data.frame(Year=Yr,
mean=xMean, data.frame(t(xCI)))
xSum
}
sumExpMeans <- sumSims(1/simMeans)
sumExpMeans4 <- sumSims(1/simMeans4)
```
These numbers are added to Figure \@ref(fig:ploty) to produce Figures \@ref(fig:plotfut) and \@ref(fig:plotfu4).
```{r plotfut, fig.cap = paste('Estimated mean time between "first tests," past and future.', NucStates)}
plotNucStatesPred <- function(x, ...){
##
## plotNucStates with future predictions
## summarized in x
##
xlim. <- range(x$Year)
ylim. <- range(
nuclearWeaponStates$yearsSinceLastFirstTest,
head(x[-1], 1), tail(x[-1], 1), na.rm=TRUE)
plotNucStates(xlim.=xlim.,
ylim.=ylim., log.='y', ...)
with(x, lines(Year, mean))
with(x, lines(Year, median, lty='dashed',
col='blue'))
with(x, lines(Year, U10, lty='dotted',
col='red'))
with(x, lines(Year, L10, lty='dotted',
col='red'))
with(x, lines(Year, U20, lty='dotted',
col='red'))
with(x, lines(Year, L20, lty='dotted',
col='red'))
legend('topleft',
c('60% and 80% confidence limits for the mean',
'mean', 'median'),
col=c('red', 'black', 'blue'),
lty=c('dotted', 'solid', 'dashed'), bty='n')
abline(h=200, lty='dotted', col='grey')
}
plotNucStatesPred(sumExpMeans)
if(save_svg){
print(getwd())
svg(paste0('NucWeaponStates_',
'FcstMeanTimeBetw1stTsts.svg'))
op <- par(cex=1.4, cex.axis=1.3,
mar=c(2,3, 2, 2)+.1)
cex2 <- 1.3
plotNucStatesPred(sumExpMeans, mtext. = FALSE,
cex.=cex2)
dev.off()
}
```
The fairly flat shape of the median and lower 10 and 20 percent lines in Figure \@ref(fig:plotfut) seem consistent with a model that is a mixture of log-normal distributions with the dominant component having a mean that is constant over time and a probability of `r 100*round(fitProlif$postprob[1], 2)` percent. The substantial curvature of the solid line forecast looks hopeful, with a mean of simulated means being almost 200 years between successive "first tests" by new nuclear-weapon states by the end of the forecasted period, `r year(tail(sumExpMeans$Year, 1))`.
The fact that the mean of the simulations exceeds the upper confidence limit for `r year(tail(sumExpMeans$Year, 1))` seems odd but can be explained by noting that this is a mixture of log-normal distributions, and the mean of a log-normal exceeds its upper quantile $q$ whenever $\sigma > 2\Phi^{-1}(q)$, where $\Phi^{-1}(q)$ = quantile $q$ of the standard normal distribution, and $\sigma$ is the standard deviation of the logarithms.^[This follows, because quantile $q$ of a log-normal is $\exp[\mu+\sigma\Phi^{-1}(q)]$ and the mean is $\exp[\mu+\sigma^2/2]$, so the mean exceeds quantile $q$ whenever $\sigma\Phi^{-1}(q) < \sigma^2/2$, i.e., when $\Phi^{-1}(q) < \sigma/2$.]
Note further that the distribution for each year in Figure \@ref(fig:plotfut) is a mixture of log-normal distributions, which means that their reciprocals, the mean numbers of "first tests" each year, will also be a mixture of log-normals with the same standard deviations on the log scale. This standard deviation is larger the farther we extrapolate into the future.
```{r plotfu4, fig.cap = paste('Estimated mean time between "first tests" considering up to a quartic model.', NucStates)}
plotNucStatesPred(sumExpMeans4)
if(save_svg){
print(getwd())
svg(paste0('NucWeaponStates_',
'QuarticFcstMeanTimeBetw1stTsts.svg'))
op <- par(cex=1.4, cex.axis=1.3,
mar=c(1,2, 0,0)+.1)
cex2 <- 1.3
plotNucStatesPred(sumExpMeans4, mtext. = FALSE,
cex.=cex2)
dev.off()
}
```
The increase over time in the _mean_ time between "first tests" in Figures \@ref(fig:plotfut) and \@ref(fig:plotfu4) suggests a desirable decrease in the rate of nuclear proliferation.
However, we are more concerned with the _shorter_ times between "first tests", and they seem all too probable, as we shall see when we simulate and `cumsum` them. To do that, we append these simulated predictions to a plot of the evolution of the number of nuclear-weapon states through the historical period.^[In these simulations, we assume a zero probability of a nuclear power giving up their nuclear weapons, even though [South Africa reportedly discontinued their nuclear weapons program in 1989](https://en.wikipedia.org/wiki/South_Africa_and_weapons_of_mass_destruction), prior to its [first universal elections in 1994](https://en.wikipedia.org/wiki/South_Africa#End_of_apartheid). We could potentially add South Africa to `nuclearWeaponStates` with the same date as Israel, then model the distribution of the time to when a nuclear-weapon state gives up its nuclear weapons using an exponential distribution. For that, we have one observed time and eight such times that are censored. Standard theory in that case says that the maximum likelihood estimate of the mean time to relinquishing nuclear weapons assuming an exponential distribution is the sum of all the times, censored or observed, divided by the number of times observed, not including the censored times in the denominator. For purposes of illustration, we will assume that South Africa dismantled its nuclear weapons 1989-12-31, though a report of an inspection by the International Atomic Energy Agency dated 1994-08-19 said they had dismantled six nuclear weapons and were still working to dismantle one more. Based on this, the mean lifetime of a nuclear-weapon state can be estimated at `r round(as.numeric(difftime(nukeDataCurrent, as.Date('1979-09-22'), units='weeks') + sum(difftime(nukeDataCurrent, nuclearWeaponStates$firstTest, units='weeks')))/52)` years. We could potentially add this to the current modeling effort, but it would not likely change the answers enough to justify the additional effort.]
```{r cumsum}
str(cumMeans <- apply(simMeans[-(1:nYrs), ],
2, cumsum))
quantile(cumMeans[nYrs,])
str(cumCI <- sumSims(
nNucStates+rbind(0, cumMeans),
pastfut$Year[-(1:(nYrs-1))]))
```
These numbers are plotted in Figure \@ref(fig:plotcum).
```{r plotcum_Fn}
plotNNucStates <- function(xfuture,
xpast=nuclearWeaponStates,
lwd.=c(1,1,1,2,2), yNPT=2.5, xlim.=NULL,
ylim.=NULL, fT_end=nukeDataCurrent, ...){
##
## plot stairsteps for xpast and lines for xfuture
## with either 5 or 7 columns in xfuture
##
if(is.null(xlim.)){
xlim. <- range(c(xpast$firstTest,
tail(xfuture$Year, 1)))
}
nColsFut <- length(xfuture)
if(is.null(ylim.))ylim. <- c(0,
tail(xfuture[[nColsFut]], 1))
plot(xlim., ylim., type='n', xlab='',
ylab='', las=1, bty='n')
##
## 1. plot xpast
##
fT_sel <- (xpast$firstTest<=fT_end)
nNucStend <- sum(fT_sel)
fT_date <- c(xpast$firstTest[1],
xpast$firstTest[fT_sel], fT_end)
# fT_date <- c(xpast$firstTest[1],
# xpast$firstTest, fT_end)
lines(fT_date,
c(0:nNucStend, nNucStend), type='s')
abline(v=NPTdate, lty='dashed', col='grey')
# abline(h=20)
xlim20. <- c(xlim.[1]-7*365, xlim.[1])
lines(xlim20., rep(20, 2), lty='dashed',
col='grey', xpd=NA)
xlim20 <- c(as.Date('1980-01-01'), xlim.[2])
lines(xlim20, rep(20, 2), lty='dashed',
col='grey')
#
text(NPTdate-.017*diff(xlim.), yNPT,
c('NPT', 'INF'), col='grey', srt=90)
##
## 2. plot xfuture
##
if(!is.null(xfuture)){
with(xfuture, lines(Year, mean),
lwd=lwd.[1])
with(xfuture, lines(Year, median,
lty='dashed', col='blue', lwd=lwd.[2]))
with(xfuture, lines(Year, U10,
lty='dotted', col='red', lwd=lwd.[3]))
with(xfuture, lines(Year, L10,
lty='dotted', col='red', lwd=lwd.[3]))
with(xfuture, lines(Year, U20,
lty='dotted', col='red', lwd=lwd.[3]))
with(xfuture, lines(Year, L20,
lty='dotted', col='red', lwd=lwd.[3]))
ncols <- 3
leg <- c(paste0('60% and 80%',
'confidence\nlimits for the mean'),
'mean', 'median')
col. <- c('red', 'black', 'blue')
#
if('predU10' %in% names(xfuture)){
leg <- c(leg, '80% prediction limits')
col. <- c(col., 'green')
with(xfuture, lines(Year, predU10,
lty='dashed', col='green',
lwd=lwd.[4]))
with(xfuture, lines(Year, predL10,
lty='dashed', col='green',
lwd=lwd.[4]))
ncols <- 4
}
if('tolU10' %in% names(xfuture)){
leg <- c(leg,
'(0.8, 0.8) tolerance limits')
col. <- c(col., 'purple')
with(xfuture, lines(Year, tolU10,
lty='dashed', col='purple',
lwd=lwd.[5]))
with(xfuture, lines(Year, tolL10,
lty='dashed', col='purple',
lwd=lwd.[5]))
ncols <- ncols+1
}
##
## 3. legend
##
lty. <- c('dotted', 'solid',
rep('dashed', 3))
# lwd. <- c(rep(par('lwd'), 3),
# rep(par('lwd')*2, 2) )
legend('topleft', leg[1:ncols],
col=col.[1:ncols], lty=lty.[1:ncols],
lwd=lwd., bty='n')
}
}
```
```{r plotcum, fig.cap = "Number of nuclear-weapon states, past and predicted mean"}
plotNNucStates(cumCI)
print(tail(cumCI, 1))
# optionally write to a file
if(save_svg){
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
png('NucWeaponStates_nucProlifPred1.png')
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1)
# const+linear
plotNNucStates(cumCI, mtext. = FALSE)
title('Number of nuclear-weapon states')
dev.off()
}
# For the 88th MORS Symposium, I want multiple subsets of this:
# 1. Up to NPT in 1970
# 2. Up to INF in 1988
# 3. Up to "Axis of Evil" in 2002
# 4. Today in 2020
# 5. Poisson data: Model as a Poisson count
# each year since 1945
# 6. BMA: 78% const.; 22% linear
# 7. To 2094 with confidence limits
if(save_svg){
# 1. NPT
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred11.png',
width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
# plotNNucStates(xfuture, xpast=nuclearWeaponStates,
# lwd.=c(1,1,1,2,2), yNPT=2.5, xlim.=NULL, ylim.=NULL, ...)
xlim. <- c(nuclearWeaponStates$firstTest[1],
tail(cumCI$Year, 1))
ylim. <- c(0, tail(cumCI$U10, 1))
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.,
fT_end=NPTdate[1])
title('Number of nuclear-weapon states')
dev.off()
}
if(save_svg){
# 2. INF
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred12.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
# plotNNucStates(xfuture, xpast=nuclearWeaponStates,
# lwd.=c(1,1,1,2,2), yNPT=2.5, xlim.=NULL, ylim.=NULL, ...)
xlim. <- c(nuclearWeaponStates$firstTest[1],
tail(cumCI$Year, 1))
ylim. <- c(0, tail(cumCI$U10, 1))
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.,
fT_end=NPTdate[2])
title('Number of nuclear-weapon states')
dev.off()
}
if(save_svg){
# 3. To "Axis of Evil"
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred13.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
# plotNNucStates(xfuture, xpast=nuclearWeaponStates,
# lwd.=c(1,1,1,2,2), yNPT=2.5, xlim.=NULL, ylim.=NULL, ...)
xlim. <- c(nuclearWeaponStates$firstTest[1],
tail(cumCI$Year, 1))
ylim. <- c(0, tail(cumCI$U10, 1))
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.,
fT_end=as.Date('2002-01-29'))
title('Number of nuclear-weapon states')
dev.off()
}
if(save_svg){
# 4. To the present
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred14.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
# plotNNucStates(xfuture, xpast=nuclearWeaponStates,
# lwd.=c(1,1,1,2,2), yNPT=2.5, xlim.=NULL, ylim.=NULL, ...)
xlim. <- c(nuclearWeaponStates$firstTest[1],
tail(cumCI$Year, 1))
ylim. <- c(0, tail(cumCI$U10, 1))
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.)
title('Number of nuclear-weapon states')
dev.off()
}
if(save_svg){
# 5. Poisson observations
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred15.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.)
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
points(nuclearWeaponStates$firstTest[-1],
rep(0.92, nNucStates-1),
col='blue', pch='.', cex=4)
tstDay. <- mean(yday(nuclearWeaponStates$firstTest))
datYrs <- (firstTstYr[1]+(1:nYrs))
tstYr <- datYrs[!(datYrs %in% firstTstYr)]
yr <- (as.Date(paste0(tstYr, '-01-01')) +
tstDay.-1)
points(yr, rep(0, nYrs-nNucStates+1),
col='blue', pch='.', cex=2)
dev.off()
}
if(save_svg){
# 5. Poisson observations with annotation
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred15a.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.)
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
points(nuclearWeaponStates$firstTest[-1],
rep(0.92, nNucStates-1),
col='blue', pch='.', cex=4)
tstDay. <- mean(yday(nuclearWeaponStates$firstTest))
datYrs <- (firstTstYr[1]+(1:nYrs))
tstYr <- datYrs[!(datYrs %in% firstTstYr)]
yr <- (as.Date(paste0(tstYr, '-01-01')) +
tstDay.-1)
points(yr, rep(0, nYrs-nNucStates+1),
col='blue', pch='.', cex=2)
text(nukeDataCurrent+100, 6,
'Model as a sequence\nof Poisson observations',
0, cex=2, col='blue')
dev.off()
}
if(save_svg){
# 6. BMA: 78% const.; 22% lin.
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred16.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.)
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
points(nuclearWeaponStates$firstTest[-1],
rep(0.92, nNucStates-1),
col='blue', pch='.', cex=4)
tstDay. <- mean(yday(nuclearWeaponStates$firstTest))
datYrs <- (firstTstYr[1]+(1:nYrs))
tstYr <- datYrs[!(datYrs %in% firstTstYr)]
yr <- (as.Date(paste0(tstYr, '-01-01')) +
tstDay.-1)
points(yr, rep(0, nYrs-nNucStates+1),
col='blue', pch='.', cex=2)
dev.off()
}
if(save_svg){
# 6. BMA: 78% const.; 22% lin. w annotations
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred16a.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
plotNNucStates(NULL, mtext. = FALSE,
lwd.=c(3,3,3,4,4), xlim. = xlim., ylim.=ylim.)
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
points(nuclearWeaponStates$firstTest[-1],
rep(0.92, nNucStates-1),
col='blue', pch='.', cex=4)
tstDay. <- mean(yday(nuclearWeaponStates$firstTest))
datYrs <- (firstTstYr[1]+(1:nYrs))
tstYr <- datYrs[!(datYrs %in% firstTstYr)]
yr <- (as.Date(paste0(tstYr, '-01-01')) +
tstDay.-1)
points(yr, rep(0, nYrs-nNucStates+1),
col='blue', pch='.', cex=2)
text(nukeDataCurrent+100, 6,
'Model as a sequence\nof Poisson observations',
0, cex=2, col='blue')
mixp <- round(100*fitProlif$postprob)
mixpt <- paste(mixp, c('const.', 'lin.'))
mixpt1 <- paste(mixpt, collapse='; ')
mixpt1b <- paste('BMA: ', mixpt1)
text(nukeDataCurrent+100, 1.5, mixpt1b,
0, cex=1.7, col='red')
dev.off()
}
if(save_svg){
# 7. Full slide
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred17.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
plotNNucStates(cumCI, mtext. = FALSE,
lwd.=c(3,3,3,4,4) )
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
yr <- as.Date(paste0(time(FirstTests$nFirstTests),
'-01-01'))
points(yr, 0.9*FirstTests$nFirstTests, col='blue',
pch='.', cex=2*(1+FirstTests$nFirstTests))
dev.off()
}
if(save_svg){
# 7. Full slide
print(getwd())
# svg('NucWeaponStates_nucProlifPred1.svg')
# svg format is not supported for Google slides
png('NucWeaponStates_nucProlifPred17a.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1, lwd=3)
# const+linear
plotNNucStates(cumCI, mtext. = FALSE,
lwd.=c(3,3,3,4,4) )
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
yr <- as.Date(paste0(time(FirstTests$nFirstTests), '-01-01'))
points(yr, 0.9*FirstTests$nFirstTests, col='blue', pch='.',
cex=2*(1+FirstTests$nFirstTests))
text(nukeDataCurrent+100, 6,
'Model as a sequence\nof Poisson observations',
0, cex=2, col='blue')
mixp <- round(100*fitProlif$postprob)
mixpt <- paste(mixp, c('const.', 'lin.'))
mixpt1 <- paste(mixpt, collapse='; ')
mixpt1b <- paste('BMA: ', mixpt1)
text(nukeDataCurrent+100, 1.5, mixpt1b,
0, cex=1.7, col='red')
dev.off()
}
```
The slopes of the mean and median lines are steeper than the recent history, but the statistical evidence does not support the naive interpretation of a slowing in nuclear proliferation that one might get from considering only the most recent data.
We repeat this analysis with the quartic BMA mixture in Figure \@ref(fig:plotcum4).
```{r plotcum4, fig.cap = "Number of nuclear-weapon states, past and predicted mean, BMA quartic fit"}
str(cumMeans4 <- apply(simMeans4[-(1:nYrs), ],
2, cumsum))
quantile(cumMeans4[nYrs,])
str(cumCI4 <- sumSims(
nNucStates+rbind(0, cumMeans4),
pastfut$Year[-(1:(nYrs-1))]))
plotNNucStates(cumCI4)
print(tail(cumCI4, 1))
# optionally write to a file
if(save_svg){
print(getwd())
svg('NucWeaponStates_nucProlifPred.svg',
height=3.5)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1,
mfrow=c(1, 2))
# const+linear
nColsFut <- length(cumCI4)
ylim. <- c(0, tail(cumCI4[[nColsFut]], 1))
plotNNucStates(cumCI, mtext. = FALSE,
ylim.=ylim.)
plotNNucStates(cumCI4, mtext. = FALSE,
ylim.=ylim.)
dev.off()
}
```
Comparing Figures \@ref(fig:plotcum) and \@ref(fig:plotcum4) shows that the higher order terms in the quartic BMA mixture widens the confidence limits, making the 10`th` percentile essentially flat with almost no additional nuclear proliferation, while the mean quickly escapes the upper limit. That sharply rising mean suggests that less than 10 percent of the simulations predict nuclear arms races that involve many nation states and many more non-state armed groups. These outcomes are not likely, but the probabilities of such outcomes seem too large to be dismissed without further consideration, especially when gambling with the future of civilization.
We replicate the simulations summarized in `cumCI4` to see how stable the numbers are for the final year in Figure \@ref(fig:plotcum4); see Figure \@ref(fig:qq4).
```{r qq4, fig.cap = "Log-normal probability plot of mean numbers of nuclear-weapon states per the quartic mixture model in the last simulated year"}
cumCI2 <- rbind(tail(cumCI, 1), tail(cumCI4, 1))
rownames(cumCI2) <- c('BMA2', 'BMA4')
for(i in 2:10){
simMeans4b <- simulate(fitProlif4, nSims, seed=i,
newdata=pastfut[2:5], type='response')
cumMeans4b <- apply(simMeans4b[-(1:nYrs), ],
2, cumsum)
cumCI4b <- sumSims(
nNucStates+rbind(0, cumMeans4b),
pastfut$Year[-(1:(nYrs-1))])
cumCI2 <- rbind(cumCI2, tail(cumCI4b, 1))
rownames(cumCI2)[i+1] <- paste0(
'BMA4', letters[i])
}
qqCI2 <- as.data.frame(
qqnorm(cumCI2$mean[-1], datax=TRUE,
log='x', ylab='simulated mean'))
with(qqCI2[1, , drop=FALSE],
points(x, y, pch=15))
with(qqCI2[1, , drop=FALSE],
text(x, y, adj=c(0, 0.5),
' <- used in other figures'))
cumCI2
```
These replications establish that the simulated `mean` number of nuclear-weapon states in the last simulated year, `r year(tail(cumCI$Year, 1))`, in Figure \@ref(fig:plotcum4) is slightly conservative relative to the simulated replicates and is definitely _not_ an outlier.
Beyond that, comparing Figures \@ref(fig:plotcum) and \@ref(fig:plotcum4) establishes that the median and lower limit are lower for the quartic BMA mixture while the upper limit is slightly higher, and the mean for the quartic BMA mixture is too large to be credible.
Ignoring the simulations of uncontrolled nuclear arms races, the median lines in Figures \@ref(fig:plotcum) and \@ref(fig:plotcum4) predict between `r round(tail(cumCI, 1)$median, 1)` and `r round(tail(cumCI4, 1)$median, 1)` at the end of the current simulated period, `r year(tail(cumCI, 1)$Year)`, adding either `r round(tail(cumCI, 1)$median-nNucStates, 1)` and `r round(tail(cumCI4, 1)$median-nNucStates, 1)` (for the constant-linear and quartic mixtures, respectively) to the current `r nNucStates` nuclear-weapon states. Those _median_ numbers are a little less than double the number of nuclear-weapon states today.
We extend this analysis by adding prediction intervals to these plots.
# Prediction limits
The simplest bounds on the future are [prediction intervals](https://en.wikipedia.org/wiki/Prediction_interval), which combine the statistical uncertainty in the estimates of mean numbers of nuclear-weapon states with the random variability in the outcomes. For this we use `rpois(., simMeanNucStByYr)` to simulate independent numbers of "first tests" by new nuclear-weapon states (or non-state entities) for each of the next `r nYrs` years in the forecast period then `cumsum` those to get simulated trajectories.
```{r simNew}
set.seed(9)
rpois. <- function(n, lambda){
##
## Some of the means are so large that
## rpois sometimes returns NAs.
## Avoid this by outputting numerics
## rather than integers in those cases.
## NOTE: This was discussed on
## r-devel@r-project.org 2020-01-19 and 20
## with the tentative conclusion that
## a change such as documented here
## might be implemented in a future version
## of R. If that happens, this
## function will no longer be needed.
##
n2 <- max(n, length(lambda))
n. <- rep_len(n, n2)
lam <- rep_len(lambda, n2)
# If Poisson mean = 0.9*.Machine$integer.max,
# an observation would have to be over
# 4600 standard deviations above the mean
# to generate an error.
big <- (lam>0.9*.Machine$integer.max)
out <- rep(NA, n2)
out[big] <- round(rnorm(sum(big),
lam[big], sqrt(lam[big])))
out[!big] <- rpois(sum(!big), lam[!big])
out
}
cumsumPred <- function(x, ...){
##
## cumsum of rpois predictions based on x
##
#
simPred <- data.frame(lapply(
x[-(1:nYrs),], rpois., n=nYrs))
cumPred <- data.frame(lapply(
simPred, cumsum))
cumPred
}
cumPred <- cumsumPred(simMeans)
cumsumC.PI <- function(cumsumPred, cumsumCI, ...){
cumPI <- sumSims(
nNucStates+rbind(0, cumsumPred),
pastfut$Year[-(1:(nYrs-1))])
prd. <- which(names(cumPI) %in%
c('L10', 'L20', 'U20', 'U10'))
names(cumPI)[prd.] <- paste0('pred',
c('L10', 'L20', 'U20', 'U10'))
# checks
dYr <- difftime(cumsumCI$Year, cumPI$Year, 'days')
if(any(as.numeric(dYr)>0))
stop('Years do not match')
rd.mean <- ((cumsumCI$mean-cumPI$mean) /
(cumsumCI$mean+cumPI$mean) )
if(any(rd.mean>0.01))
stop('means do not match')
# cbind
cumC.PI <- cbind(cumsumCI, cumPI[prd.])
cumC.PI
}
cumC.PI <- cumsumC.PI(cumPred, cumCI)
```
We add this to the image in Figure \@ref(fig:plotcum) to create Figure \@ref(fig:plotPred).
```{r plotPred, fig.cap = "Number of nuclear-weapon states, past and predicted"}
plotNNucStates(xfuture=cumC.PI)
```
We do the same for the quartic BMA model in Figure \@ref(fig:plotcum4) to produce Figure \@ref(fig:plotPred4). In both Figures \@ref(fig:plotPred) and \@ref(fig:plotPred4) the most likely scenarios, especially the median line and the space between the 60 percent confidence limits, predict a continuation of nuclear proliferation. It's difficult to imagine how that could continue without also substantively increasing the risk of nuclear war and therefore also of the extinction of civilization.
```{r plotPred4, fig.cap = "Number of nuclear-weapon states, past and predicted per quartic BMA model"}
cumPred4 <- cumsumPred(simMeans4)
cumC.PI4 <- cumsumC.PI(cumPred4, cumCI4)
plotNNucStates(xfuture=cumC.PI4)
# optionally write to a file
if(save_svg){
print(getwd())
svg('NucWeaponStates_nucProlifPredInt.svg',
height=3.5)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1,
mfrow=c(1, 2))
# const+linear
nColsFut <- length(cumC.PI4)
ylim. <- c(0, tail(cumC.PI4[[nColsFut]], 1))
plotNNucStates(cumC.PI, mtext. = FALSE,
ylim.=ylim.)
plotNNucStates(cumC.PI4, mtext. = FALSE,
ylim.=ylim.)
dev.off()
}
```
We can also summarize the simulations in `cumPred` to estimate the probabilities of having 1, 2, 3, 4, and 5 new nuclear weapon states for each year in the prediction period between `r simStart` and `r simEnd` in Figure \@ref(fig:prolif). This is another way of evaluating the sensibility of pretending there will be no further nuclear proliferation: Not likely.
```{r prolif, fig.cap = "Probabilities of the time to the next 5 new nuclear-weapon states using the constant-linear mixture model"}
plotProbs <- function(x, label_year=2050,
cex.txt = 1, ...){
# adj.=matrix(c(.53, -.4, .55, 1.15), 2, byrow=TRUE), ...){
##
## Probability distribution of the next 1:5
## new nuclear-weapon states
##
## from x = cumPred, a data.frame
##
maxNewNucSt <- 5
probs <- function(x, n=maxNewNucSt){
p <- colMeans(outer(x, 0:(n-1), '>'))
p
}
probProlif. <- apply(as.matrix(x), 1, probs)
probProlif <- ts(t(probProlif.),
currentYear+1)
colnames(probProlif) <- 1:maxNewNucSt
ylims <- probProlif.[c(1, maxNewNucSt), c(1, nYrs)]
midLine <- round(mean(ylims[c(1, 4)]), 2)
futYrs <- time(probProlif)
plot(range(futYrs), range(probProlif), type='n',
xlab='', las=1, ylab='', ...)
# matplot(time(probProlif), probProlif,
# type='l', xlab='', las=1, ylab='', ...)
for(iNNS in 1:maxNewNucSt){
lines(as.numeric(futYrs), probProlif[, iNNS],
lty=iNNS, col=iNNS)
}
abline(h=midLine, lty='dotted', ...)
#
midTime <- rep(NA, maxNewNucSt)
for(i in 1:maxNewNucSt){
it <- min(which(probProlif[,i]>=midLine))
# cat(time(probProlif)[it], '')
midTime[i] <- time(probProlif)[it]
text(midTime[i], midLine,
paste0(midTime[i], '\n', i),
cex=cex.txt)
}
text(currentYear+0.95*nYrs, midLine,
paste0('p =\n', midLine), cex=cex.txt)
# lab3l_year
abline(v=label_year, lty='dotted')
iyr <- which(time(probProlif) == label_year)
label_p <- probProlif[iyr, ]
lines(time(probProlif)[c(1, iyr)],
rep(label_p[1], 2), lty='dotted', ...)
text(time(probProlif)[1+iyr/4], label_p[1],
paste0('p =\n', round(label_p[1], 2)), cex=cex.txt)
#
list(midLine=midLine, midTime=midTime,
probProlif=probProlif,
label_year=c(label_year=label_year,
label_prob=label_p))
}
probProlif <- plotProbs(cumPred)
# optionally write to a file
if(save_svg){
print(getwd())
# svg('NucWeaponStates_nucProbs.svg',
# height=3.5)
png('NucWeaponStates_nucProbs1.png', width=960)
op <- par(cex=1.8, cex.axis=1.8, lwd=3,
mar=c(2,3, 2, 0)+.1)
# const+linear
cex.txt <- 1.6
ylim. <- c(0, tail(
probProlif$probProlif[, 1], 1))
plotProbs(cumPred, ylim=ylim.,
cex.txt=cex.txt)
dev.off()
}
tail(probProlif$probProlif, 1)
```
Figure \@ref(fig:prolif) says there is a `r round(100*tail(probProlif$probProlif[, 1], 1))` percent chance of another nuclear-weapon state by `r tail(time(probProlif$probProlif), 1)` with a `r 100*probProlif$midLine` percent chance of at least 1 by `r probProlif$midTime[1]`.
We replicate Figure \@ref(fig:prolif) for the quartic BMA model in Figure \@ref(fig:prolif4).
```{r prolif4, fig.cap = "Probabilities of the time to the next 5 new nuclear-weapon states using the quartic BMA model"}
probProlif4 <- plotProbs(cumPred4)
tail(probProlif4$probProlif, 1)
# optionally write to a file
if(save_svg){
print(getwd())
svg('NucWeaponStates_nucProbs.svg',
height=3.5)
op <- par(cex=1, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1,
mfrow=c(1, 2))
# const+linear
cex.txt <- 0.7
ylim. <- c(0, tail(
probProlif$probProlif[, 1], 1))
plotProbs(cumPred, ylim=ylim.,
cex.txt=cex.txt)
plotProbs(cumPred4, ylim=ylim.,
cex.txt=cex.txt)
dev.off()
}
```
The quartic BMA model summarized in Figure \@ref(fig:prolif4) is more optimistic than the constant-linear BMA model summarized in Figure \@ref(fig:prolif): The estimated probability of one more nuclear-weapon state by 2050 drops from `r round(window(probProlif$probProlif, 2050, 2050)[,1], 2)` for the constant-linear mixture to `r round(window(probProlif4$probProlif, 2050, 2050)[,1], 2)` for the quartic BMA model.
The conclusions from both models include the following:
\begin{center}
\textbf{The current structure of international relations}
\textbf{seems to threaten the extinction of civilization}
\end{center}
To better quantify the uncertainty in modeling, we can also construct tolerance intervals for the time to the next new nuclear-weapon state.
# Tolerance limits
We want to add tolerance limits to Figures \@ref(fig:plotPred) and \@ref(fig:plotPred4) in addition to the prediction limits. [A (p, 1-$\alpha$) tolerance interval provides limits within which at least a proportion (p) of the population falls with a given level of confidence (1−$\alpha$).](https://en.wikipedia.org/wiki/Tolerance_interval)
We compute this as follows:
1. Find the L10 and U90, the central 80 percent confidence limits for the Poisson means for individual years summarized in Figures \@ref(fig:plotcum) and \@ref(fig:plotcum4).
2. Simulate Poisson distributed random variables about each of those limits for each individual year in the future.
3. Then we `cumsum` those simulations.
4. Select TL10 and TU90 to include at least 80 percent of `cumsum`med Poisson simulations about each of L10 and U90 individually. We do this as follows:
4.1. Set the initial value for `TU90` as the 80th percentile of the distributions of the simulated trajectories around `U90`. Similarly set the initial value for `TL10` as the 20`th` percentile of the distributions of the simulated trajectories around `L10`.
4.2. Compute the coverage probabilities of (TL10`, `TU90`) of the simulations around `L10` and `U90`, respectively. If both those coverage probabilities exceed 80 percent, quit.
4.3. Else, compute the coverage probabilities of (`max(0, TL10-1), TU90`) and (`TL10, TU90+1`) under both `L10` and `U90`, respectively. If both of those coverage probabilities exceed 80 percent, quit with the smallest of the two. If only one exceeds 80 percent, quit with that one. If neither exceeds 80 percent, replace (`TL10, TU90`) with the max of the two [or with (TL10, TU90+) if they are both the same] and repeat this analysis.
```{r simTol}
cumsumTol <- function(x=cumCI, ...){
##
## cumsum of rpois predictions based on x
##
# simTolU. <- rpois.(nSims*nYrs, diff(x$U10))
simTolU <- matrix(rpois.(nSims*nYrs,
diff(x$U10)), nYrs)
simTolL <- matrix(rpois.(nSims*nYrs,
diff(x$L10)), nYrs)
#
cumTolU <- apply(simTolU, 2, cumsum)
cumTolL <- apply(simTolL, 2, cumsum)
#
cumTolU10 <- apply(cumTolU, 1,
quantile, probs=0.8)
cumTolL10 <- apply(cumTolL, 1,
quantile, probs=0.2)
# check coverage for each year
# with both limits in case too many
# upper sims are below cumTolL10 and
# too many lower sims are above
# cumTolU10
In <- ((cumTolL10 <= cumTolL) &
(cumTolU <= cumTolU10))
pIn <- apply(In, 1, mean)
pOut <- (pIn < 0.8)
while(any(pOut)){
stop('Problem in tolSim: ',
'pIn < 0.8; need algorithm to fix;',
' not debugged, because it',
' seems not to have actually',
' happened.\n')
# The following not debugged:
cumTolU10.1 <- (cumTolU10[pOut]+1)
cumTolL10.1 <- pmax(0,
cumTolL10[pOut]-1)
# Need In.1 and In.9
In.1 <- ((cumTolL10.1 <=
cumTolL[pOut, drop=FALSE])
& (cumTolU[pOut,, drop=FALSE]
<= cumTolU10[pOut]))
In.9 <- ((cumTolL10[pOut] <=
cumTolL[pOut, drop=FALSE])
& (cumTolU[pOut,, drop=FALSE]
<= cumTolU10.1))
pIn.1 <- apply(In.1, 1, mean)
pIn.9 <- apply(In.9, 1, mean)
pO.1 <- (pIn.1 < 0.8)
pO.9 <- (pIn.9 < 0.8)
# If both pIn.1 and pIn.9 >=0.8,
# pick the smallest
# else if one is, use that,
# else iterate
# ... the following is not correct,
# partly because neither cumTolL10
# nor cumTolU10 have been changed
pOut[pOut] <- pO
}
#
x$tolL10 <- (nNucStates + c(0, cumTolL10))
x$tolU10 <- (nNucStates + c(0, cumTolU10))
#
as.data.frame(x)
}
cumTol <- cumsumTol(cumC.PI)
```
The results for the constant-linear mixture appear in Figure \@ref(fig:plotTol).
```{r plotTol, fig.cap = "Number of nuclear-weapon states with prediction and tolerance limits"}
plotNNucStates(cumTol)
if(save_svg){
print(getwd())
# svg('NucWeaponStates_nucProlifTolInt.svg',height=3.5)
png('NucWeaponStates_nucProlifTolInt1.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 1)+.1, lwd=3)
# const+linear
# ylimTol1 <- c(0, tail(cumCI$U10, 1))
plotNNucStates(cumTol, mtext. = FALSE,
lwd.=c(3,3,3,4,4) )
# lwd.=c(3,3,3,4,4), ylim.=ylimTol1)
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
yr <- as.Date(paste0(time(FirstTests$nFirstTests), '-01-01'))
points(yr, 0.9*FirstTests$nFirstTests, col='blue', pch='.',
cex=2*(1+FirstTests$nFirstTests))
dev.off()
}
if(save_svg){
print(getwd())
# svg('NucWeaponStates_nucProlifTolInt.svg',height=3.5)
png('NucWeaponStates_nucProlifTolInt1a.png', width=960)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 1)+.1, lwd=3)
# const+linear
# ylimTol1 <- c(0, tail(cumCI$U10, 1))
plotNNucStates(cumTol, mtext. = FALSE,
lwd.=c(3,3,3,4,4) )
# lwd.=c(3,3,3,4,4), ylim.=ylimTol1)
title('Number of nuclear-weapon states')
# lines(c(nuclearWeaponStates$firstTest[1], Today),
# c(0,0), lwd=3, col='red')
# points(nuclearWeaponStates$firstTest[-1], 0,
# nuclearWeaponStates$firstTest[-1], 1,
# cex=3, col='red')
yr <- as.Date(paste0(time(FirstTests$nFirstTests), '-01-01'))
points(yr, 0.9*FirstTests$nFirstTests, col='blue', pch='.',
cex=2*(1+FirstTests$nFirstTests))
text(nukeDataCurrent+100, 6,
'Model as a sequence\nof Poisson observations',
0, cex=2, col='blue')
text(nukeDataCurrent+100, 1.5,
'BMA: 79% const.; 21% lin.',
0, cex=1.7, col='red')
dev.off()
}
```
The upper limit line in Figure \@ref(fig:plotTol) is higher than that in Figure \@ref(fig:plotcum). It gives us a bit more humility regarding the value of current knowledge. However, the difference is not enough to substantively alter our conclusions, namely that nuclear proliferation is likely and should not be ignored.
Do we get the same considering the quartic BMA model? See Figure \@ref(fig:plotTol4).
```{r plotTol4, fig.cap = "Number of nuclear-weapon states with prediction and tolerance limits per quartic BMA model"}
cumTol4 <- cumsumTol(cumC.PI4)
plotNNucStates(cumTol4)
if(save_svg){
print(getwd())
svg('NucWeaponStates_nucProlifTolInt.svg',
height=3.5)
op <- par(cex=1.7, cex.axis=1.4,
mar=c(2,3, 2, 0)+.1,
mfrow=c(1, 2))
# const+linear
nColsFut <- length(cumTol4)
ylim. <- c(0, tail(cumTol4[[nColsFut]], 1))
plotNNucStates(cumTol, mtext. = FALSE,
ylim.=ylim.)
plotNNucStates(cumTol4, mtext. = FALSE,
ylim.=ylim.)
dev.off()
}
```
Indeed, the conclusion from Figure \@ref(fig:plotTol4) is the same as before: Nuclear proliferation is likely to continue until something makes it impossible for anyone to make more nuclear weapons for a very long time.
<!--This might happen either as a result of (a) a nuclear war leading to the destruction of civilization or (b) a major and unprecedented political movement that strengthens international law to the point that the poor, weak, and disfranchised have effective judicial redress for grievances.-->
# Discussion
A growing number of leading figures have said that as long as the world maintains large nuclear arsenals, it is only a matter of time before there is a nuclear war. Concerns like this have been expressed by two former US Secretaries of Defense ([Robert McNamara](https://en.wikipedia.org/wiki/Robert_McNamara)^[@McNamara:2003] and [William Perry](https://en.wikipedia.org/wiki/William_Perry)), two former US Secretaries of State ([Henry Kissinger](https://en.wikipedia.org/wiki/Henry_Kissinger) and [George Schultz](https://en.wikipedia.org/wiki/George_Shultz)), former US Senator [Sam Nunn](https://en.wikipedia.org/wiki/Sam_Nunn),^[@Shultz:2019] and others with, for example, the [Nuclear Threat Initiative](https://en.wikipedia.org/wiki/Nuclear_Threat_Initiative). [Daniel Ellsberg](https://en.wikipedia.org/wiki/Daniel_Ellsberg) has said that a nuclear war will most likely generate a nuclear winter that lasts several years during which 98 percent of humanity will starve to death if they do not die of something else sooner.^[@Ellsberg:2017b]
Banerjee and Duflo, two of the three who won the 2019 Nobel Memorial Prize in Economics, have noted that neither economic nor political stability are assured for any country, including the United States, China and India. In particular, they predict that economic growth will almost certainly slow substantially in the latter two, leaving many poor people in desperate economic straits.^[@Banerjee:2019. Various journalists and academic researchers have expressed concern about increases in ethic violence in various countries and whether electoral transitions of power will continue, even in the US. See, e.g., @Klaas:2019] Internal problems in the US, China, India or any other nuclear-weapon state could push political leaders to pursue increasingly risky foreign adventures, like Argentina did in 1982,^[@Wikip:Faulklands] possibly leading to a war that could produce nuclear Armageddon.[^WikivTime2]
The evidence compiled in the present work only seems to increase the urgency of limiting the threat of nuclear war and nuclear proliferation in particular.
In the 20 years following the first test of a nuclear weapon on `r nuclearWeaponStates$firstTest[1]` by the US, four more nations acquired such weapons. In the 50 years since the Non-Proliferation Treaty took effect in 1970, another four acquired them.^[This uses a commonly accepted list of existing nuclear-weapon states and when they each first tested a nuclear weapon. The sources for this are contained in the help file for the [`nuclearWeaponStates`](https://www.rdocumentation.org/packages/Ecdat/versions/0.3-7/topics/nuclearWeaponStates) dataset in the [Ecfun](https://www.rdocumentation.org/packages/Ecdat/versions/0.3-7) package for R.] Our analysis of the available data considering only the dates of these first tests suggests that nuclear proliferation may have been slowing throughout this period. However, that apparent trend was not statistically significant in the model we fit.
Bayesian Model Averages (BMA) is known to generally produce better predictions than single model fits. Accordingly, we've estimated confidence, prediction, and tolerance limits for the number of new nuclear-weapon states `r nYrs` years into the future based on two BMA models with mixtures of either a constant with a linear model or a constant with terms up to quartic in the time since the very first test of a nuclear weapon.
We can expect that some non-nuclear nations and terrorist groups would eagerly pursue nuclear weapons if such seemed feasible unless some unprecedented change in international law provided them with effective nonviolent recourse to perceived threats.
Moreover, these weapons will likely become more available with the passage of time unless (a) a nuclear war destroys everyone's ability to make more such weapons for a long time, or (b) an international movement has far more success than similar previous efforts in giving the poor, weak and disfranchised effective nonviolent means for pursuing a redress of grievances.
# Acknowledgments
The authors wish to thank Val Nereo for suggestions that seem to have improved this article. As usual, the authors retain responsibility for any remaining errors and poorly presented observations.
# References
[^quickNukes]: In addition to the 32 currently non-nuclear-weapon states with "sufficient fissile material to make nuclear weapons if they wished", per @Toon:2007, the inspector general of the US Department of Energy concluded in 2009 (in its most recent public accounting) that enough highly enriched uranium was missing from US inventories to make at least five nuclear bombs comparable to those that destroyed substantial portions of Hiroshima and Nagasaki in 1945. Substantially more weapons-grade material may be missing in other countries, especially Russia (@Malone:2018).
[^WikivTime]: For precursors to the current study that involve censored estimation of time to a nuclear war, see @Wikiv:extinct and @Wikiv:Armageddon.
[^WikivTime2]: The risks of a nuclear Armageddon are documented in a series of simulations published in refereed academic journals, each more detailed and more disconcerting than the previous. All assume that many firestorms will be produced, because (a) the areas targeted will likely be much more susceptible to firestorms than the underground or isolated sites used to test nuclear weapons, and (b) many of the weapons used will have yields substantially greater than those employed in Hiroshima and Nagasaki. For a discussion of that literature, see @Wikiv:extinct and @Wikiv:Armageddon.
[^KP]: @Levy:2007 say that Gallucci was a special adviser on WMDs to US Presidents Clinton and G. W. Bush. The Wikipedia article on him says he was US Assistant Secretary of State for Political-Military Affairs from July 13, 1992 to October 11, 1994 under Presidents George H. W. Bush and Bill Clinton but not G. W. Bush. Later, per @Gallucci:2001, "In March 1998, the Department of State announced his appointment as Special Envoy to deal with the threat posed by the proliferation of ballistic missiles and weapons of mass destruction. He held this position until January 2001." G. W. Bush became US President 2001-01-20. Thus, if Gallucci served under G. W. Bush, it was only for a few days. Similar remarks about the US helping Pakistan's nuclear program were made by [Richard Barlow](https://en.wikipedia.org/wiki/Richard_Barlow_(Intelligence_analyst)#cite_note-wp-1), a CIA analyst who reported these questionable activities to a committee of the US House as noted by @Levy:2007. Barlow was reportedly severely punished for honestly answering questions in a classified briefing to an oversight committee of the US House. Barlow said that US assistance to Pakistan's nuclear weapons program was in exchange for Pakistan's help in supplying rebels in Afghanistan fighting Soviet occupation. This was during the [Iran-Contra affair](https://en.wikipedia.org/wiki/Iran%E2%80%93Contra_affair), which exposed actions of officials of the Reagan administration to pursue foreign policy objectives in Central America in blatant violation of law passed by Congress and signed by the President.
[^Pak2]: @Burr:2012, @Burr:2013. There have also been reports that China helped Pakistan obtain nuclear weapons. However, [China has vigorously denied those charges, many if not all of which may not be credible, having originated with the US government.](https://en.wikipedia.org/wiki/Pakistan_and_weapons_of_mass_destruction) See @Wikip:PakNuc.
[^Confidence_Interval]:["Confidence intervals"](https://en.wikipedia.org/wiki/Confidence_interval) bound the predicted mean number of nuclear-weapon states for each future year considered. Central 80 percent ["prediction intervals"](https://en.wikipedia.org/wiki/Prediction_interval) are limits that include the central 80 percent of distribution of the number of nuclear-weapon states. They add the uncertainty in the modeled Poisson process to the uncertainty of estimating the mean of that process for each future year considered. We will also compute
(0.8, 0.8) ["tolerance intervals"](https://en.wikipedia.org/wiki/Tolerance_interval#Relation_to_other_intervals); $(p, 1-\alpha)$ tolerance intervals have a probability of $(1-\alpha)$ of containing a proportion of at least $p$ of all future observations.
nuc-references.bib
[edit | edit source]@online{Aftergood:2006,
title = {DPRK: Nuclear Weapons Program},
url = {https://fas.org/nuke/guide/dprk/nuke/},
urldate = {2020-02-05},
publisher = {Federation of American Scientists}
}
@Book{Bacevich:2008,
author={Bacevich, Andrew},
title = {The Limits of Power: The End of American Exceptionalism},
publisher={Metropolitan Books},
year = {2008}
}
@Book{Banerjee:2019,
author={Abhijit V. Banerjee and Esther Duflo},
title = {Good Economics for Hard Times},
publisher={Public Affairs},
year = {2019}
}
@Article{Bazzi:2019,
author = {Bazzi, Mohamad},
title = {Both Saudi Arabia and the United States Are Probably Guilty of War Crimes in Yemen},
date = {2019-05-17},
url = {https://www.thenation.com/article/war-crimes-united-states-saudi-arabia-yemen/},
urldate = {2019-12-23},
journal = {The Nation}
}
@Book{Benjamin:2016,
author={Benjamin, Medea},
title={Kingdom of the Unjust: Behind the US-Saudi Connection},
publisher={ORBooks},
year={2016},
url={https://en.wikipedia.org/wiki/Kingdom_of_the_Unjust},
urldate = {2019-05-25}
}
@Book{Berton:1980,
author = {Berton, Pierre},
title = {Invasion of Canada, 1812-1813},
publisher={McClelland and Stewart},
year = {1980}
}
@Article{Borger:2010,
author = {Borger, Julian},
title = {Pakistan nuclear weapons at risk of theft by terrorists, US study warns},
journal = {The Guardian},
date = {2010-04-12},
url = {https://www.theguardian.com/world/2010/apr/12/pakistan-nuclear-weapons-security-fears},
urldate = {2019-05-01}
}
@Book{BoxDraper:1987,
author = {George E. P. Box and Norman Draper},
title = {Empirical Model-Building and Response Surfaces},
publisher = {Wiley},
year = {1987}
}
@Book{Burnham:1998,
author = {Kenneth P. Burnham and David R. Anderson},
title = {Model Selection and Multimodel Inference: A practical information-theoretic approach},
publisher={Springer},
year = {1998}
}
@Book{Burnham:2002,
author = {Kenneth P. Burnham and David R. Anderson},
title = {Model Selection and Multimodel Inference: A practical information-theoretic approach},
publisher={Springer},
year = {2002},
related={Model Selection and Multimodel Inference: A practical information-theoretic approach},
relatedstring={See :},
}
@Article{Burr:2012,
author = {Burr, William},
title = {New Documents Spotlight Reagan-era Tensions over Pakistani Nuclear Program},
date = {2012-04-25},
url = {https://www.wilsoncenter.org/publication/new-documents-spotlight-reagan-era-tensions-over-pakistani-nuclear-program},
publisher={Wilson Center},
urldate = {2019-05-01}
}
@Article{Burr:2013,
author = {Burr, William},
title = {Pakistan's Illegal Nuclear Procurement Exposed in 1987: Arrest of Arshed Pervez Sparked Reagan Administration Debate over Sanctions},
date = {2013-11-22},
url = {https://nsarchive2.gwu.edu/nukevault/ebb446/},
publisher={National Security Archives},
urldate = {2019-05-01}
}
@online{Bush:2002,
author = {George W. Bush},
title={President Delivers State of the Union Address},
url={https://georgewbush-whitehouse.archives.gov/news/releases/2002/01/print/20020129-11.html},
date={2002-01-29},
urldate={2019-12-23},
publisher={Whitehouse Archives}
}
@Book{Claeskens:2008,
author = {Gerda Claeskens and Nils Lid Hjort},
title = {Model selection and model averaging},
publisher={Cambridge U. Pr.},
year = {2008}
}
@Book{Chomsky:2017,
author = {Noam Chomsky},
title = {Who rules the world?},
publisher={Picador},
year = {2017}
}
@Article{Cohen:2009,
author = {Cohen, Stephen P.},
title = {The US-Pakistan Strategic Relationship and Nuclear
Safety/Security},
date = {2008-06-12},
url = {https://www.brookings.edu/testimonies/the-u-s-pakistan-strategic-relationship-and-nuclear-safetysecurity},
publisher={Brookings Institution},
urldate = {2019-05-01}
}
@online{CRS:2016,
title = {Comprehensive Nuclear-Test-Ban Treaty: Background and Current Developments},
date = {2016-09-01},
url = {https://www.everycrsreport.com/files/20160901_RL33548_9db23ae2fefe501277b8b3bfc505870c36aced66.html#_Toc461031566},
publisher = {Congressional Research Service, Library of Congress, Washington, DC},
urldate = {2020-02-02}
}
@Book{Ellsberg:2017,
author = {Ellsberg, Daniel},
title = {Doomsday Machine: Confessions of a nuclear war planner},
publisher={Bloomsbury},
year = {2017},
url = {https://en.wikipedia.org/wiki/Daniel_Ellsberg#The_Doomsday_Machine},
urldate = {2019-05-15}
}
@Article{Ellsberg:2017b,
author = {Daniel Ellsberg and Amy Goodman and Juan González},
title = {Daniel Ellsberg Reveals He was a Nuclear War Planner, Warns of Nuclear Winter & Global Starvation},
date = {2017-12-06},
journal = {Democracy Now!},
url = {https://www.democracynow.org/2017/12/6/doomsday_machine_daniel_ellsberg_reveals_he},
urldate = {2019-12-23}
}
@Book{Fogelsong:1995,
author={Fogelsong, David S.},
title = {America's Secret War against Bolshevism: {U.S.} Intervention in the Russian Civil War, 1917-1920},
publisher={UNC Press},
year = {1995}
}
@Article{Forbes:2006,
author = {Forbes},
title = {AFX News Limited: "Saudia Arabia working on secret nuclear program with Pakistan help - report"},
journal={Forbes},
date={2006-03-28},
url={https://web.archive.org/web/20120115022055/http://www.forbes.com/feeds/afx/2006/03/28/afx2629000.html},
urldate={2019-11-29}
}
@Article{Fromkin:2006,
author = {David Fromkin},
title = {Stuck in the Canal},
journal = {New York Times},
date = {2006-10-28},
url = {https://www.nytimes.com/2006/10/28/opinion/28fromkin.html},
urldate = {2020-02-03}
}
@Article{Gallucci:2001,
author = {Robert Gallucci},
title = {Robert L. Gallucci: Dean, School of Foreign Service, Georgetown University},
date = {2001-10-15},
url = {http://globetrotter.berkeley.edu/people2/Gallucci/gallucci-vita.html},
urldate = {2019-05-06}
}
@online{GlobalSec:Uriso,
title = {Weapons of Mass Destruction ({WMD}): Uranium Isotopes},
url = {https://www.globalsecurity.org/wmd/intro/u-isotopes.htm},
urldate={2020-02-05},
publisher = {GlobalSecurity.org}
}
@techreport{Graham:2003,
author = {Bob Graham and Porter Goss and Richard Shelby and Nancy Pelosi},
title = {Joint inquiry into intelligence community activities before and after the terrorist attacks of September 11, 2001: Report of the {U.S.} Senate Select Committee on Intelligence and {U.S.} House Permanent Select Committee on Intelligence, December 2002 with letter of transmittal to the Director of National Intelligence},
Date = {2003-01-29},
publisher = {Intelligence Committee of the {U.S.} House of Representatives},
address = {Washington, DC, USA},
url = {https://web.archive.org/web/20160715183528/http://intelligence.house.gov/sites/intelligence.house.gov/files/documents/declasspart4.pdf},
urldate = {2019-12-23}
}
@techreport{Halperin:1966,
author = {Morton Halperin},
title = {The 1958 Taiwan Straits Crisis: A documentary history},
Date = {1966-12},
publisher = {RAND Corporation},
url = {https://www.rand.org/content/dam/rand/pubs/research_memoranda/2006/RM4900.pdf},
urldate = {2019-12-24}
}
@Article{Helfand:2019,
author = {Ira Helfand},
title = {Ban the Bomb–Before Our Luck Runs Out},
date = {2019-06-01},
url = {https://progressive.org/magazine/ban-the-bomb-helfand/},
urldate = {2019-11-03},
journal = {The Progressive}
}
@Vignette{Helske:2017,
author = {Helske, Jouni and Matti Vihola},
title = {bssm: Bayesian Inference of Non-linear and
Non-Gaussian State Space Models in R},
date = {2017-01-04},
url = {https://www.rdocumentation.org/packages/bssm/versions/0.1.7/vignettes/bssm.Rmd},
urldate= {2019-09-01}
}
@Article{Klaas:2019,
author = {Brian Klaas},
title = {Everyone knows the 2020 election will be divisive. But will it also be violent?},
journal = {Washington Post},
date = {2019-09-05},
url = {https://www.washingtonpost.com/opinions/2019/09/05/everyone-knows-election-will-be-divisive-will-it-also-be-violent/},
urldate = {2019-12-23}
}
@incollection{Kolko:1968,
author = {Kolko, Gabriel},
title = {Report on the destruction of dikes: Holland 1944-45 and Korea 1953},
editor = {John Duffett},
booktitle = {Against the crime of silence: Proceedings of the Russell International War Crimes Tribunal},
publisher = {O'Hare},
address = {New York},
year = {1968},
pages = {224-226}
}
@Article{Levy:2007,
author = {Levy,Adrian and Cathy Scott-Clark},
title = {Pakistan nuclear weapons at risk of theft by terrorists, US study warns},
journal = {The Guardian},
date = {2007-10-13},
url = {https://www.theguardian.com/world/2007/oct/13/usa.pakistan},
urldate = {2019-05-01}
}
@Article{Malone:2018,
author = {Malone, Patrick and Jeffrey Smith},
title = {Plutonium is missing, but the government says nothing: Losses of civilian nuclear material are usually disclosed but when the government loses nuclear bomb ingredients it stays mum},
journal = {Center for Public Integrity},
date = {2018-07-16},
url = {https://publicintegrity.org/national-security/plutonium-is-missing-but-the-government-says-nothing/},
urldate = {2019-05-01}
}
@Book{McNamara:2003,
author = {Robert McNamara and James G. Blight},
title = {Wilson's ghost: Reducing the risk of conflict, killing, and catastrophe in the 21st century},
publisher={Public Affairs},
year = {2003}
}
@Article{MEM:SaudiNuc,
author={{Middle East Monitor}},
title={Israel ‘is selling nuclear information’ to Saudi Arabia},
journal={Middle East Monitor},
date={2018-05-31},
url={https://www.middleeastmonitor.com/20180531-israel-is-selling-nuclear-information-to-saudi-arabia/},
urldate={2019-11-29}
}
@Article{OConnor:2019,
author = {O'Connor, Tom},
title = "Turkey Has {U.S.} Nuclear Weapons, Now It Says It Should Be Allowed to Have Some of Its Own",
journal = {Newsweek},
date = {2019-09-04},
url = {https://www.newsweek.com/turkey-us-nuclear-weapons-its-own-1457734},
urldate = {2019-10-17}
}
@online{ORB:2008,
author = {{Opinion Research Business}},
title={Update on Iraqi Casualty Data},
url={https://web.archive.org/web/20080201144355/http://www.opinion.co.uk/Newsroom_details.aspx?NewsId=88},
date = {2018-01},
urldate={2019-12-23},
publisher={Opinion Research Business}
}
@online{osti:espionage,
author = {OSTI},
title={Espionage and the Manhattan Project (1940-1945)},
publisher={{Office of Scientific and Technical Information}, US Department of Energy},
url = {https://www.osti.gov/opennet/manhattan-project-history/Events/1942-1945/espionage.htm},
urldate={2020-01-08}
}
@Article{Pierson:2017,
author = {Charles Pierson},
title = {The Atomic Bomb and the First Korean War},
journal = {CounterPunch},
date = {2017-09-08},
url = {https://www.counterpunch.org/2017/09/08/the-atomic-bomb-and-the-first-korean-war/},
urldate = {2019-12-24}
}
@incollection{Raftery:1995,
author = {Raftery, Adrian},
title = {Bayesian model selection in social research (with Discussion)},
editor = {Peter V. Marsden},
booktitle = {Sociological Methodology 1995},
publisher = {Blackwells},
address = {Cambridge, Mass.},
year = {1995},
pages = {111-196}
}
@Book{Rhodes:1986,
author={Richard Rhodes},
title = {The Making of the Atomic Bomb},
publisher={Simon & Schuster},
year = {1986},
url = {https://en.wikipedia.org/wiki/The_Making_of_the_Atomic_Bomb},
urldate = {2019-12-23}
}
@online{Riedel:2008,
author = {Bruce Riedel},
title = {Saudi Arabia: Nervously Watching Pakistan},
publisher = {Brookings},
date = {2008-01-28},
url = {https://web.archive.org/web/20120205150154/http://www.brookings.edu/opinions/2008/0128_saudi_arabia_riedel.aspx},
urldate = {2020-02-02}
}
@Article{Riedel:2012,
author = {Bruce Riedel},
title = {JFK's Overshadowed Crisis},
journal = {The National Interest},
date = {2012-06-28},
url = {https://nationalinterest.org/article/jfks-overshadowed-crisis-7073},
urldate={2020-02-04}
}
@Article{Shultz:2019,
author = {Shultz, George P. and William J. Perry and Sam Nunn},
title = {The Threat of Nuclear War Is Still With Us},
date = {2019-04-10},
url = {https://www.wsj.com/articles/the-threat-of-nuclear-war-is-still-with-us-11554936842},
urldate = {2019-11-03},
journal = {Wall Street Journal}
}
@Article{Sood:2008,
author = {Sood, Vikram},
title = {America fails the IQ test},
date = {2008-01-17},
url = {https://www.hindustantimes.com/india/america-fails-the-iq-test/story-G0wNJq151e1Rs9fzATbmEN.html},
urldate = {2020-04-20},
journal = {Hindustani Times}
}
@Article{Spokane:1976,
author={Spokane Daily Chronicle},
title = {Student Designs Nuclear Bomb},
journal={Spokane Daily Chronicle},
date = {1976-10-09},
url = {https://news.google.com/newspapers?id=_cESAAAAIBAJ&sjid=4fgDAAAAIBAJ&pg=2166,2187720&dq=john-aristotle-phillips&hl=en},
urldate={2020-01-08}
}
@Article{Toksabay:2019,
author = {Toksabay, Ece},
title = {Erdogan says it's unacceptable that Turkey can't have nuclear weapons},
journal = {Reuters},
date = {2019-09-04},
url = {https://www.reuters.com/article/us-turkey-nuclear-erdogan-idUSKCN1VP2QN},
urldate = {2019-10-17}
}
@Article{Toon:2017,
author = {Owen B. Toon and Alan Robock and Michael Mills and Lili Xia},
title = {Asia treds the nuclear path, unaware that self-assured destruction would result from nuclear war},
journal = {Journal of Asian Studies},
date = {May 2017},
url = {http://climate.envsci.rutgers.edu/pdf/ToonAsianStudies.pdf},
urldate={2020-02-29},
pages = {437-456},
volume = {76},
number = {2}
}
@Article{Toon:2007,
author = {Toon, O. B. and R. P. Turco and A. Robock and C. Bardeen and L. Oman and G. L. Stenchikov},
title = {Atmospheric effects and societal consequences of regional scale nuclear conflicts and acts of individual nuclear terrorism},
journal = {Atmospheric Chemistry and Physics},
pages = {1973-2002},
volume = {7},
number = {8},
date = {2007-04-19},
url = {http://climate.envsci.rutgers.edu/pdf/acp-7-1973-2007.pdf},
urldate = {2019-10-17}
}
@online{UN:1970,
author={{United Nations Office for Disarmament Affairs}},
title = {Treaty on the Non-Proliferation of Nuclear Weapons},
date={1970-03-05},
url = {http://disarmament.un.org/treaties/t/npt},
urldate={2020-01-08},
publisher={United Nations Office for Disarmament Affairs}
}
@online{Wikip:28pages,
title={The 28 pages},
url={https://en.wikipedia.org/wiki/The_28_pages#External_links},
urldate={2019-12-23},
publisher={Wikipedia}
}
@online{Wikip:AlliesInRussia,
title={Allied intervention in the Russian Civil War},
url={https://en.wikipedia.org/wiki/Allied_intervention_in_the_Russian_Civil_War},
urldate={2019-12-23},
publisher={Wikipedia}
}
@online{Wikip:ArabIsrael,
title = {Arab-Israeli conflict},
url = {https://en.wikipedia.org/wiki/Arab%E2%80%93Israeli_conflict},
urldate={2020-02-07},
publisher={Wikipedia}
}
@online{Wikip:AxisOfEvil,
title={Axis of evil},
url={https://en.wikipedia.org/wiki/Axis_of_evil},
urldate={2019-12-23},
publisher={Wikipedia}
}
@online{Wikip:Bhabha,
title = {Homi J. Bhabha},
url = {https://en.wikipedia.org/wiki/Homi_J._Bhabha#Death},
urldate={2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:ChinaNuc,
title={China and weapons of mass destruction},
url={https://en.wikipedia.org/wiki/China_and_weapons_of_mass_destruction},
urldate={2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:Ellsberg,
title={Daniel Ellsberg},
url={https://en.wikipedia.org/wiki/Daniel_Ellsberg},
urldate={2019-12-24},
publisher={Wikipedia}
}
@online{Wikip:EnrichedU,
title = {Enriched Uranium},
url = {https://en.wikipedia.org/wiki/Enriched_uranium#Highly_enriched_uranium_(HEU)},
urldate={2020-02-05},
publisher={Wikipedia}
}
@online{Wikip:Faulklands,
title={Falklands War},
url={https://en.wikipedia.org/wiki/Falklands_War},
urldate={2019-12-23},
publisher={Wikipedia}
}
@online{Wikip:FranceNuc,
title = {France and weapons of mass destruction},
url = {https://en.wikipedia.org/wiki/France_and_weapons_of_mass_destruction},
urldate = {2020-04-21},
publisher = {Wikipedia}
}
@online{Wikip:HistNucWeapons,
title={History of nuclear weapons},
url={https://en.wikipedia.org/wiki/History_of_nuclear_weapons},
urldate={2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:IndiaNuc,
title={India and weapons of mass destruction},
url={https://en.wikipedia.org/wiki/India_and_weapons_of_mass_destruction},
urldate={2020-02-04},
publisher = {Wikipedia}
}
@online{Wikip:IraqWar,
title={Casualties of the Iraq War},
url={https://en.wikipedia.org/wiki/Casualties_of_the_Iraq_War},
urldate={2019-12-23},
publisher={Wikipedia}
}
@online{Wikip:IsraelNuc,
title={Nuclear weapons and Israel},
url = {https://en.wikipedia.org/wiki/Nuclear_weapons_and_Israel},
urldate={2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:Manhattan,
title={Manhattan Project},
url = {https://en.wikipedia.org/wiki/Manhattan_Project},
urldate={2020-02-05},
publisher = {Wikipedia}
}
@online{Wikip:ModelsWrong,
title={All models are wrong},
url={https://en.wikipedia.org/wiki/All_models_are_wrong},
urldate = {2019-12-22},
publisher={Wikipedia}
}
@online{Wikip:NPT,
title={Treaty on the Non-Proliferation of Nuclear Weapons},
url={https://en.wikipedia.org/wiki/Treaty_on_the_Non-Proliferation_of_Nuclear_Weapons},
urldate={2020-01-08},
publisher={Wikipedia}
}
@online{Wikip:NKoreaNuc,
title = {North Korea and weapons of mass destruction},
url = {https://en.wikipedia.org/wiki/North_Korea_and_weapons_of_mass_destruction},
urldate={2020-02-05},
publisher={Wikipedia}
}
@online{Wikip:NKoreaNuc06,
title={2006 North Korean nuclear test},
url={https://en.wikipedia.org/wiki/2006_North_Korean_nuclear_test},
urldate={2020-02-02},
publisher={Wikipedia}
}
@online{Wikip:NucPhys,
title={Nuclear physics},
url = {https://en.wikipedia.org/wiki/Nuclear_physics},
urldate = {2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:NucTimeline,
title={Timeline of nuclear weapons development},
url = {https://en.wikipedia.org/wiki/Timeline_of_nuclear_weapons_development},
urldate={2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:PakNuc,
title={Pakistan and weapons of mass destruction},
url={https://en.wikipedia.org/wiki/Pakistan_and_weapons_of_mass_destruction},
urldate = {2020-01-08},
publisher={Wikipedia}
}
@online{Wikip:SaudiNuc,
title={Nuclear program of Saudi Arabia},
url={https://en.wikipedia.org/wiki/Nuclear_program_of_Saudi_Arabia},
urldate={2019-05-26},
publisher={Wikipedia}
}
@online{Wikip:SaudiYemen,
title={Saudi Arabian-led intervention in Yemen},
url={https://en.wikipedia.org/wiki/Saudi_Arabian-led_intervention_in_Yemen},
urldate = {2019-12-23},
publisher={Wikipedia}
}
@online{wikip:Soviet,
title={Soviet atomic bomb project},
url={https://en.wikipedia.org/wiki/Soviet_atomic_bomb_project},
urldate={2020-01-08},
publisher={Wikipedia}
}
@online{Wikip:studentBomb,
title={John Aristotle Phillips},
url={https://en.wikipedia.org/wiki/John_Aristotle_Phillips},
urldate={2020-02-08}
}
@online{Wikip:TaiwanStrait,
title={First Taiwan Strait Crisis},
url={https://en.wikipedia.org/wiki/First_Taiwan_Strait_Crisis#Aftermath:_China_and_nuclear_weapons},
urldate={2020-02-04},
publisher={Wikipedia}
}
@online{Wikip:Trinity,
title={Trinity (nuclear test)},
url={https://en.wikipedia.org/wiki/Trinity_(nuclear_test)},
urldate = {2019-12-23},
publisher={Wikipedia}
}
@online{Wikip:Ures,
title = {List of countries by uranium reserves},
url = {https://en.wikipedia.org/wiki/List_of_countries_by_uranium_reserves},
urldate={2020-02-05},
publisher={Wikipedia}
}
@online{Wikip:Wiener,
title={Wikipedia, "Wiener process"},
url={https://en.wikipedia.org/wiki/Wiener_process},
urldate = {2019-12-23},
publisher={Wikipedia}
}
@online{Wikiv:WinTerror,
title={Winning the war on terror},
url={https://en.wikiversity.org/wiki/Winning_the_War_on_Terror},
urldate = {2019-05-25},
publisher={Wikiversity}
}
@online{Wikiv:extinct,
title={Time to the extinction of civilization},
url={https://en.wikiversity.org/wiki/Time_to_extinction_of_civilization},
urldate = {2019-11-29},
publisher={Wikiversity}
}
@online{Wikiv:Armageddon,
title={Time to nuclear Armageddon},
url={https://en.wikiversity.org/wiki/Time_to_nuclear_Armageddon},
urldate = {2019-05-25},
publisher={Wikiversity}
}
@online{Wolfram:Wiener,
title={Wolfram, "Wiener Process"},
url={http://mathworld.wolfram.com/WienerProcess.html},
urldate = {2019-12-23},
publisher={Wolfram}
}
@Article{Yusuf:2016,
author = {Yusuf, Moeed},
title = {An India-Pakistan Crisis: Should we care?},
date = {2016-11-29},
url = {https://warontherocks.com/2016/11/an-indian-pakistan-crisis-should-we-care/},
urldate = {2016-11-29},
journal = {War on the Rocks}
}