Bubble script for /u/swolebird

### We need the following packages for this example.
packages <- c("lubridate", "fredr", "shape", "forecast",
              "ggthemes", "zoo", "tidyverse", "tis", "showtext", "tsibble",
              "tsibbledata", "feasts", "fable", "dplyr", "tsbox",
              "here", "magick", "magrittr", "reshape2", "gtrendsR",
              "scales")

### Install packages if needed, then load them quietly
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
invisible(lapply(packages, "library", quietly = TRUE,
                 character.only = TRUE, warn.conflicts = FALSE))

### Load API Keys
api_key_fred <- "INSERT_YOUR_FRED_API_KEY_HERE"

# I like the "Economist" ggplot theme, so use it by default.
# theme_set(theme_economist() + scale_colour_economist())
theme_set(theme_economist())

####################################################################
### Add recession bars to ggplot graphs
####################################################################
geom_recession_bars <- function (date_start, date_end) {

  date_start = as.Date(date_start, origin="1970-01-01")
  date_end   = as.Date(date_end,   origin="1970-01-01")

  recessions_tibble <- tibble(

    peak = as.Date(
      c("1857-06-01", "1860-10-01", "1865-04-01", "1869-06-01",
        "1873-10-01", "1882-03-01", "1887-03-01", "1890-07-01",
        "1893-01-01", "1895-12-01", "1899-06-01", "1902-09-01",
        "1907-05-01", "1910-01-01", "1913-01-01", "1918-08-01",
        "1920-01-01", "1923-05-01", "1926-10-01", "1929-08-01",
        "1937-05-01", "1945-02-01", "1948-11-01", "1953-07-01",
        "1957-08-01", "1960-04-01", "1969-12-01", "1973-11-01",
        "1980-01-01", "1981-07-01", "1990-07-01", "2001-03-01",
        "2007-12-01")),

    trough = as.Date(
      c("1858-12-01", "1861-06-01", "1867-12-01", "1870-12-01",
        "1879-03-01", "1885-05-01", "1888-04-01", "1891-05-01",
        "1894-06-01", "1897-06-01", "1900-12-01", "1904-08-01",
        "1908-06-01", "1912-01-01", "1914-12-01", "1919-03-01",
        "1921-07-01", "1924-07-01", "1927-11-01", "1933-03-01",
        "1938-06-01", "1945-10-01", "1949-10-01", "1954-05-01",
        "1958-04-01", "1961-02-01", "1970-11-01", "1975-03-01",
        "1980-07-01", "1982-11-01", "1991-03-01", "2001-11-01",
        "2009-06-01")
    )
  )

  recessions_trim <- recessions_tibble %>%
    filter(peak   >= min(date_start) &
             trough <= max(date_end))

  if (nrow(recessions_trim) > 0) {

    recession_bars = geom_rect(data        = recessions_trim,
                               inherit.aes = F,
                               fill        = "darkgray",
                               alpha       = 0.25,
                               aes(xmin = as.Date(peak,   origin="1970-01-01"),
                                   xmax = as.Date(trough, origin="1970-01-01"),
                                   ymin = -Inf, ymax = +Inf))
  } else {

    recession_bars = geom_blank()
  }

}

### Set my FRED API key to access the FRED database.
### The actual value is in Include/API_Keys.R
### You may request an API key at:
### https://research.stlouisfed.org/useraccount/apikeys
fredr_set_key(api_key_fred)

### What date ranges do we want?
date_start <- "1952-01-01"
date_end   <- "2019-04-01"

### A frequency selector to automate selecting the data frequency
#frequency <- c("m", 12, "Seasonal12") # For monthly data
frequency <- c("q",  4, "Seasonal4")  # For quarterly data
#frequency <- c("sa", 2, "Seasonal2")  # For semiannual data
#frequency <- c("a",  1, "Seasonal1")  # For annual data

###########################################################################
### Import Household Wealth from FRED
###########################################################################
series <- "TNWBSHNO"

data <- fredr(series_id = series, frequency = frequency[1],
              observation_start = as.Date(date_start),
              observation_end   = as.Date(date_end)) %>%
        as_tsibble(index = "date")

date             <- data %>% pull("date")
household_wealth <- data %>% pull("value")
household_wealth_decomposed <- household_wealth %>%
                               ts(frequency = as.numeric(frequency[2])) %>%
                               mstl()

household_wealth_seasonal  <- household_wealth_decomposed %>% seasonal()
household_wealth_trend     <- household_wealth_decomposed %>% trendcycle()
household_wealth_remainder <- household_wealth_decomposed %>% remainder()

### You can use trend, trend + remainder, or trend + seasonal + remainder
household_wealth <- household_wealth_trend + household_wealth_remainder

###########################################################################
### Import Nominal GDP from FRED
###########################################################################
series <- "GDP"

data <- fredr(series_id = series, frequency = frequency[1],
              observation_start = as.Date(date_start),
              observation_end   = as.Date(date_end)) %>%
        as_tsibble(index = "date")

nominal_gdp <- data %>% pull("value")
nominal_gdp_decomposed <- nominal_gdp %>%
                          ts(frequency = as.numeric(frequency[2])) %>%
                          mstl()

nominal_gdp_seasonal  <- nominal_gdp_decomposed %>% seasonal()
nominal_gdp_trend     <- nominal_gdp_decomposed %>% trendcycle()
nominal_gdp_remainder <- nominal_gdp_decomposed %>% remainder()

### You can use trend, trend + remainder, or trend + seasonal + remainder
nominal_gdp  <- nominal_gdp_trend + nominal_gdp_remainder

###########################################################################
### Import Real GDP from FRED
###########################################################################
series <- "GDPC1"

data <- fredr(series_id = series, frequency = frequency[1],
              observation_start = as.Date(date_start),
              observation_end   = as.Date(date_end)) %>%
        as_tsibble(index = "date")

real_gdp <- data %>% pull("value")
real_gdp_decomposed <- real_gdp %>%
                       ts(frequency = as.numeric(frequency[2])) %>%
                       mstl()

real_gdp_seasonal  <- real_gdp_decomposed %>% seasonal()
real_gdp_trend     <- real_gdp_decomposed %>% trendcycle()
real_gdp_remainder <- real_gdp_decomposed %>% remainder()

### You can use trend, trend + remainder, or trend + seasonal + remainder
real_gdp  <- real_gdp_trend + real_gdp_remainder

### Graph globals
annotationcolor <- "black"

### Common caption strings to make doing the caption easier
real_gdp_cap <- "U.S. Bureau of Economic Analysis, Real Gross Domestic Product [GDPC1]\n"
nom_gdp_cap  <- "U.S. Bureau of Economic Analysis, Gross Domestic Product [GDP]\n"
wealth_cap   <- "Board of Governors of the Federal Reserve System (US), Households and nonprofit organizations; net worth, Level [TNWBSHNO]\n"

nw_cap <- paste(wealth_cap,
                nom_gdp_cap,
                "retrieved from FRED, Federal Reserve Bank of St. Louis")

nwr_cap <- paste(wealth_cap,
                 nom_gdp_cap,
                 real_gdp_cap,
                "retrieved from FRED, Federal Reserve Bank of St. Louis")

##########################################################################
### Graph 1: Household Wealth graphed alongside Nominal GDP
##########################################################################

### Normalize both Household Wealth and Nominal GDP with respect to
### the start date (1952-01-01)
norm_household_wealth <- household_wealth / household_wealth[1]
norm_nominal_gdp      <- nominal_gdp      / nominal_gdp[1]

### Create a dataframe for graphing
data_df <- data.frame(date, norm_household_wealth, norm_nominal_gdp)

### Create the graph
title    <- paste("Household Wealth and Nominal GDP growth relative to",
                  date_start)
subtitle <- "Recessions marked with vertical bars"
caption  <- nw_cap
xlab     <- "Year"
ylab     <- paste("Growth relative to ", date_start)

### Arrow settings for the annotations
segment_df_1 <- data.frame(x1 = as.Date("1995-01-01"),
                           x2 = as.Date("1998-01-01"), y1 = 32, y2 = 32)
segment_df_2 <- data.frame(x1 = as.Date("2000-07-01"),
                           x2 = as.Date("2003-07-01"), y1 = 52, y2 = 52)
segment_df_3 <- data.frame(x1 = as.Date("2012-01-01"),
                           x2 = as.Date("2015-01-01"), y1 = 82, y2 = 82)

p <- ggplot(data = data_df,
            mapping = aes(x = date, y = norm_nominal_gdp)) +

   theme_economist() + scale_colour_economist() +
   theme(legend.title = element_blank()) +
   labs(title    = title,
        subtitle = subtitle,
        caption  = caption,
        x        = xlab,
        y        = ylab) +

   geom_line(data = data_df,
             size = 1.3,
             aes(y     = norm_nominal_gdp,
                 color = "Nominal GDP")) +

   geom_line(data = data_df,
             size = 1.3,
             aes(y     = norm_household_wealth,
                 color = "Household Wealth")) +

   # Add recession bars
   geom_recession_bars(min(date), max(date)) +

   annotate("text",
            x     = as.Date("1990-01-01"),
            y     = 32,
            label = "Dot.Com\nStock Bubble",
            size  = 4,
            color = annotationcolor) +
   geom_segment(data  = segment_df_1,
                aes(x = x1, y = y1, xend = x2, yend = y2),
                arrow = arrow(length = unit(0.03, "npc")),
                size  = 1.3,
                color = annotationcolor) +

   annotate("text",
            x     = as.Date("1997-01-01"),
            y     = 52,
            label = "Housing\nBubble",
            size  = 4,
            color = annotationcolor) +
   geom_segment(data = segment_df_2,
                aes(x = x1, y = y1, xend = x2, yend = y2),
                arrow = arrow(length = unit(0.03, "npc")),
                size  = 1.3,
                color = annotationcolor) +

   annotate("text",
            x     = as.Date("2009-01-01"),
            y     = 82,
            label = "Current\nBubble",
            size  = 4,
            color = annotationcolor) +
   geom_segment(data = segment_df_3, aes(x = x1, y = y1, xend = x2, yend = y2),
                arrow = arrow(length = unit(0.03, "npc")),
                size  = 1.3,
                color = annotationcolor)
print(p)

readline(prompt = "Press [ENTER} to continue...")

###########################################################################
### Graph 2: Household Wealth to GDP ratio
###########################################################################

### Compute the ratio of Household Wealth to Nominal GDP
wealth_to_gdp    <- household_wealth   / nominal_gdp

### Create a dataframe for graphing
data_df <- data.frame(date, wealth_to_gdp)

### Add lines to show median (1951-1995) +/- 1 stdev
sd      <- wealth_to_gdp[date <= 1995] %>% sd()
mean    <- wealth_to_gdp[date <= 1995] %>% mean()
mean_p1 <- mean + sd
mean_m1 <- mean - sd

### If you plot wealth_to_gdp, you see there are some long-term trends
### that sit beneath the bubbles.  We're primarily interest in the bubbles
### as they tend not to fall below those trends when they pop.  I model
### the trend as three piecewise sections, take the local minima inside
### each trend, and use a linear regression to determine the floors.
### Then I use the "segments()" command to plot them.
###
### Pre-1960:         Wealth_Floor <- -72.60170 + 0.038941*Year
### 1960-1978.75:     Wealth_Floor <- +52.16769 - 0.024760*Year
### 1982 -> Present:  Wealth_Floor <- -63.48539 + 0.033682*Year

### Create the graph
title    <- "Ratio of Household Wealth to Nominal DDP"
subtitle <- "Blue lines show 1952-1995 mean +/- 1 std dev\nRecessions marked with vertical bars"
caption  <- nw_cap
xlab     <- "Year"
ylab     <- "Ratio"

### Arrow settings for the annotations
segment_df_1 <- data.frame(x1 = as.Date("1995-01-01"),
                           x2 = as.Date("1998-01-01"), y1 = 4.3, y2 = 4.3)
segment_df_2 <- data.frame(x1 = as.Date("2000-07-01"),
                           x2 = as.Date("2003-07-01"), y1 = 4.8, y2 = 4.8)
segment_df_3 <- data.frame(x1 = as.Date("2012-01-01"),
                           x2 = as.Date("2015-01-01"), y1 = 5.2, y2 = 5.2)

### Trend line for 1952 -1958, 1958 - 1978, and 1978 - present
segment_df_trend_1 <- data.frame(x1 = as.Date("1952-01-01"),
                                 x2 = as.Date("1958-09-03"),
                                 y1 = 3.411132, y2 = 3.67096)
segment_df_trend_2 <- data.frame(x1 = as.Date("1958-09-03"),
                                 x2 = as.Date("1978-12-14"),
                                 y1 = 3.670960, y2 = 3.16885)
segment_df_trend_3 <- data.frame(x1 = as.Date("1978-12-14"),
                                 x2 = as.Date(now()),
                                 y1 = 3.168855, y2 = 4.52617)

### Generate the plot
p <- ggplot(data = data_df,
            mapping = aes(x = date, y = wealth_to_gdp)) +

   theme_economist() + scale_colour_economist() +
   theme(legend.position = "none") +
   labs(title    = title,
        subtitle = subtitle,
        caption  = caption,
        x        = xlab,
        y       = ylab) +

   geom_line(data = data_df,
             size = 1.3,
             aes(y = wealth_to_gdp, color = "Wealth to GDP")) +

   # Add recession bars
   geom_recession_bars(min(date), max(date)) +

   geom_hline(yintercept = mean,
              size = 1,
              linetype = "solid",
              alpha = 1,
              color = "red") +

   geom_hline(yintercept = mean_p1,
              size = 1,
              linetype = "dotted",
              alpha = 1,
              color = "red") +

   geom_hline(yintercept = mean_m1,
              size = 1,
              linetype = "dotted",
              alpha = 1,
              color = "red") +

   # Draw the three trend lines
   geom_segment(data = segment_df_trend_1,
                size = 1.3,
                aes(x = x1, y = y1, xend = x2, yend = y2, color = "Segments")) +
   geom_segment(data = segment_df_trend_2,
                size = 1.3,
                aes(x = x1, y = y1, xend = x2, yend = y2, color = "Segments")) +
   geom_segment(data = segment_df_trend_3,
                size = 1.3,
                aes(x = x1, y = y1, xend = x2, yend = y2, color = "Segments")) +

   # Annotate the graphs
   annotate("text", x = as.Date("1990-01-01"), y = 4.3,
            label = "Dot.Com\nStock Bubble",
            size  = 5, color = annotationcolor) +
   geom_segment(data = segment_df_1, arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2), size = 1.3,
                color = annotationcolor) +

   annotate("text", x = as.Date("1997-01-01"), y = 4.8,
            label = "Housing\nBubble",
            size  = 5, color = annotationcolor) +
   geom_segment(data = segment_df_2, arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2), size = 1.3,
                color = annotationcolor) +

   annotate("text", x = as.Date("2009-01-01"), y = 5.2,
            label = "Current\nBubble",
            size  = 5, color = annotationcolor) +
   geom_segment(data = segment_df_3, arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2), size = 1.3,
                color = annotationcolor)

print(p)

readline(prompt = "Press [ENTER} to continue...")

###########################################################################
### Graph 3: Bubble Size to GDP
###########################################################################

### Subtract out the long-term trend mentioned in the previous section
### to get the size of the asset bubbles compared to GDP

asset_bubbles <- rep(NA, length(wealth_to_gdp))
early <- date <= as.Date("1959-09-03")
late  <- date >  as.Date("1978-12-14")
mid   <- !early & !late
asset_bubbles[early] <- wealth_to_gdp[early] +
                              72.6017 - 0.038941 * decimal_date(date[early])

asset_bubbles[mid]   <- wealth_to_gdp[mid]   -
                              52.1677 + 0.024760 * decimal_date(date[mid])

asset_bubbles[late]  <- wealth_to_gdp[late]  +
                              63.4854 - 0.033682 * decimal_date(date[late])

### Add lines to show the mean and +/1 one std.dev line for 1952-1995
mean    <- asset_bubbles[date <= 1995] %>% mean()
sd      <- asset_bubbles[date <= 1995] %>% sd()
mean_p1 <- mean + sd
mean_m1 <- mean - sd

data_df <- data.frame(date, asset_bubbles)

### Graph start
title    <- "Ratio of Asset Bubble Size to GDP"
subtitle <- "Blue lines show 1952-1995 mean +/- 1 std dev\nRecessions marked with vertical bars"
caption  <- nw_cap
xlab     <- "Year"
ylab     <- "Multiple of GDP"

### Arrows for the annotations
segment_df_1 <- data.frame(x1 = as.Date("1995-01-01"),
                           x2 = as.Date("1998-01-01"), y1 = 0.45, y2 = 0.45)
segment_df_2 <- data.frame(x1 = as.Date("2000-07-01"),
                           x2 = as.Date("2003-07-01"), y1 = 0.75, y2 = 0.75)
segment_df_3 <- data.frame(x1 = as.Date("2012-01-01"),
                           x2 = as.Date("2015-01-01"), y1 = 0.75, y2 = 0.75)

### Draw the graph
p <- ggplot(data = data_df, mapping = aes(x = date, y = asset_bubbles)) +
   theme_economist() + scale_colour_economist() +
   theme(legend.position = "none") +
   geom_line(data = data_df,
             aes(y = asset_bubbles, color = "Asset Bubbles"), size = 1.3) +
  # Add recession bars
  geom_recession_bars(min(date), max(date)) +

   geom_hline(yintercept = mean,    size = 1, linetype = "solid",
              alpha = 0.5, color = "red") +
   geom_hline(yintercept = mean_p1, size = 1, linetype = "solid",
              alpha = 0.2, color = "red") +
   geom_hline(yintercept = mean_m1, size = 1, linetype = "solid",
              alpha = 0.2, color = "red") +
   labs(title = title, subtitle = subtitle, caption = caption,
        x = xlab, y = ylab) +
   annotate("text", x = as.Date("1990-01-01"), y = 0.45,
            label = "Dot.Com\nStock Bubble",
            size  = 5, color = annotationcolor) +
   geom_segment(data = segment_df_1, arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2), size = 1.3,
                color = annotationcolor) +

   annotate("text", x = as.Date("1997-01-01"), y = 0.75,
            label = "Housing\nBubble",
            size  = 5,  color = annotationcolor) +
   geom_segment(data = segment_df_2, arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2), size = 1.3,
                color = annotationcolor) +

   annotate("text", x = as.Date("2009-01-01"), y = 0.75,
            label = "Current\nBubble",
            size  = 5, color = annotationcolor) +
   geom_segment(data = segment_df_3, arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2), size = 1.3,
                color = annotationcolor)
print(p)

readline(prompt = "Press [ENTER} to continue...")

###########################################################################
### Graph 4: Bubble Size in real 2019 US dollars
###########################################################################

### FRED Real GDP gives Real GDP in 2012 dollars.  To get it in real
### 2019 dollars, we need to multiple by 1.12 to account for inflation
inflation <- 1.12

### Calulate the real 2019 dollar value of the asset bubbles
asset_bubbles_real_dollars <- inflation * asset_bubbles * real_gdp / 1000

### Our data frame for graphing...
data_df <- data.frame(date, asset_bubbles_real_dollars)

### Create the graph
title    <- "Asset Bubble Size in Real 2019 US$"
subtitle <- "Recessions marked with vertical bars"
caption  <- nwr_cap
xlab     <- "Year"
ylab     <- "Trillions of 2019 US $"

### Arrows for the annotations
segment_df_1 <- data.frame(x1 = as.Date("1995-01-01"),
                           x2 = as.Date("1998-01-01"), y1 = 7, y2 = 7)
segment_df_2 <- data.frame(x1 = as.Date("2001-01-01"),
                           x2 = as.Date("2004-01-01"), y1 = 13, y2 = 13)
segment_df_3 <- data.frame(x1 = as.Date("2013-01-01"),
                           x2 = as.Date("2016-01-01"), y1 = 16, y2 = 16)

### Draw the graph
p <- ggplot(data = data_df,
            mapping = aes(x = date, y = asset_bubbles_real_dollars)) +

  # Graph settings
  theme_economist() +
  scale_colour_economist() +
  theme(legend.position = "none") +
  labs(title    = title,
       subtitle = subtitle,
       caption  = caption,
       x        = xlab,
       y        = ylab) +

  geom_line(data  = data_df,
            size  = 1.3,
            aes(y     = asset_bubbles_real_dollars,
                color = "Asset Bubbles")) +

  # Add recession bars
  geom_recession_bars(min(date), max(date)) +

   annotate("text",
            x     = as.Date("1991-01-01"),
            y     = 7,
            label = "Dot.Com\nStock Bubble",
            size  = 5,
            color = annotationcolor) +

   geom_segment(data  = segment_df_1,
                arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2),
                size  = 1.3,
                color = annotationcolor) +

   annotate("text",
            x     = as.Date("1998-01-01"),
            y     = 13,
            label = "Housing\nBubble",
            size  = 5,
            color = annotationcolor) +

   geom_segment(data  = segment_df_2,
                arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2),
                size  = 1.3,
                color = annotationcolor) +

   annotate("text",
            x     = as.Date("2010-01-01"),
            y     = 16,
            label = "Current\nBubble",
            size  = 5,
            color = annotationcolor) +

   geom_segment(data  = segment_df_3,
                arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2),
                size  = 1.3,
                color = annotationcolor)
print(p)

readline(prompt = "Press [ENTER} to continue...")

###########################################################################
### Graph 5: Bubble Size in real 2019 US dollars (zoom on 1990 - present)
###########################################################################

### This graph is just a zoomed-in version of the graph in the previous
### section, so all we have to do is add a "xlim()" to the appropriate
### range.

segment_df_1 <- data.frame(x1 = as.Date("1996-06-01"),
                           x2 = as.Date("1998-01-01"), y1 = 7, y2 = 7)
segment_df_2 <- data.frame(x1 = as.Date("2002-01-01"),
                           x2 = as.Date("2004-01-01"), y1 = 13, y2 = 13)
segment_df_3 <- data.frame(x1 = as.Date("2014-01-01"),
                           x2 = as.Date("2016-01-01"), y1 = 16, y2 = 16)

p <- ggplot(data = data_df,
            mapping = aes(x = date,
                          y = asset_bubbles_real_dollars)) +

   # Graph settings
   theme_economist() +
   scale_colour_economist() +
   theme(legend.position = "none") +
   labs(title    = title,
        subtitle = subtitle,
        caption  = caption,
        x        = xlab,
        y        = ylab) +

   geom_line(data = data_df,
             size = 1.3,
             aes(y = asset_bubbles_real_dollars,
                 color = "Asset Bubbles")) +

  # Add recession bars
  geom_recession_bars(min(date), max(date)) +


   scale_x_date(date_breaks = "2 years",
                labels = date_format("%Y"),
                limits = as.Date(c("1990-01-01", now()))) +

   annotate("text",
            x     = as.Date("1994-01-01"),
            y     = 7,
            label = "Dot.Com\nStock Bubble",
            size  = 5,
            color = annotationcolor) +

   geom_segment(data  = segment_df_1,
                arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2),
                size  = 1.3,
                color = annotationcolor) +

   annotate("text",
            x     = as.Date("2000-01-01"),
            y     = 13,
            label = "Housing\nBubble",
            size  = 5,
            color = annotationcolor) +

   geom_segment(data  = segment_df_2,
                arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2),
                size  = 1.3,
                color = annotationcolor) +
   annotate("text",
            x     = as.Date("2012-01-01"),
            y     = 16,
            label = "Current\nBubble",
            size  = 5,
            color = annotationcolor) +

   geom_segment(data  = segment_df_3,
                arrow = arrow(length = unit(0.03, "npc")),
                aes(x = x1, y = y1, xend = x2, yend = y2),
                size  = 1.3,
                color = annotationcolor)

print(p)