library(quantmod)


install.packages("psymonitor")

knitr::opts_chunk$set(eval = TRUE, echo = TRUE, warning = FALSE,
                      message = FALSE, comment = NA)

# CRAN Packages 
# -------------
# If any of these packages are not installed, write
# install.packages("<name of package>")
# e.g., install.packages("MultipleBubbles")

library(psymonitor)  # For testting for bubble monitoring
library(ggplot2)     # To handle plots
library(knitr)       # for nice looking tables

# For the illustration purposes we will use data on the 
# credit risk in the European sovereign sector, that is 
# proxied by an index constructed as a GDP weighted 
# 10-year government bond yield of the GIIPS (Greece, 
# Ireland, Italy, Portugal, and Spain) countries, and 
# comes with the 'psymonitor' package.

# European credit risk

library(psymonitor)
data(spread)
head(spread)

ggplot(spread, aes(date, value)) +
  geom_line() +
  ggtitle("Figure 1: Credit Risk in the European Sovereign Sector") +
  xlab("Year") +
  ylab("Index")

# Next, define a few parameters for the test and the simulation.

y        <- spread$value
obs      <- length(y)
swindow0 <- floor(obs * (0.01 + 1.8 / sqrt(obs))) # set minimal window size
IC       <- 2  # use BIC to select the number of lags
adflag   <- 6  # set the maximum nuber of lags to 6
yr       <- 2  
Tb       <- 12*yr + swindow0 - 1  # Set the control sample size
nboot    <- 99  # set the number of replications for the bootstrap

# Next, estimate the PSY test statistic using PSY() and its corresponding 
# bootstrap-based critical values using cvPSYwmboot().

bsadf <- PSY(y, swindow0 = swindow0, IC = IC,
         adflag = adflag)  # estimate the PSY test statistics sequence

quantilesBsadf <- cvPSYwmboot(y, swindow0 = swindow0, IC = IC,
                  adflag = adflag, Tb = Tb, nboot = 99,
                  nCores = 2) # simulate critical values via wild bootstrap. Note that the number of cores is arbitrarily set to 2.

# Next, identify crisis periods, defined as periods where the test 
# statistic is above its corresponding critical value, using the locate()
# function.

dim          <- obs - swindow0 + 1 
monitorDates <- spread$date[swindow0:obs]
quantile95   <- quantilesBsadf %*% matrix(1, nrow = 1, ncol = dim)
ind95        <- (bsadf > t(quantile95[2, ])) * 1
periods      <- locate(ind95, monitorDates)  # Locate crisis periods

# Finally, print a table that holds the identified crisis periods 
# with the help of the disp() function.

crisisDates <- disp(periods, obs)  #generate table that holds crisis periods

print(crisisDates)

ggplot() + 
  geom_rect(data = crisisDates, aes(xmin = start, xmax = end, 
                                    ymin = -Inf, ymax = Inf), alpha = 0.5) + 
  geom_line(data = spread, aes(date, value)) +
  labs(title = "Figure 2: Crisis episodes in the European Sovereign Sector",
       subtitle = "June 1997 - June 2016",
       caption = "Notes: The solid line is the 10-year government bond yield spread between the PIIGS countries and Germany and the shaded areas are the periods where the PSY statistic exceeds its 95% bootstrap critical value.", 
       x = "Year", y = "Index") 

###

# SP 500 dividend yield

#library(psymonitor)  # For testting for bubble monitoring
#library(ggplot2)     # To handle plots
#library(knitr)       # for nice looking tables

data(snp)
snp$pd <-  1/snp$value
head(snp)

y        <- snp$pd
obs      <- length(y)
r0       <- 0.01 + 1.8/sqrt(obs)
swindow0 <- floor(r0*obs)
dim      <- obs - swindow0 + 1
IC       <- 2
adflag   <- 6
yr       <- 2
Tb       <- 12*yr + swindow0 - 1
nboot    <- 99

# Next, we run the test and simulate critical 
# values via the bootstrap.

bsadf <- PSY(y, swindow0, IC, adflag)
quantilesBsadf <- cvPSYwmboot(y, swindow0, IC, adflag, Tb, nboot, nCores = 2) #Note that the number of cores

monitorDates <- snp$date[swindow0:obs]
quantile95   <- quantilesBsadf %*% matrix(1, nrow = 1, ncol = dim)
ind95        <- (bsadf > t(quantile95[2, ])) * 1
periods      <- locate(ind95, monitorDates)

bubbleDates <- disp(periods, obs)
kable(bubbleDates, caption = "Bubble and Crisis Periods in the S&P 500")

ggplot() + 
  geom_rect(data = bubbleDates, aes(xmin = start, xmax = end, 
                                    ymin = -Inf, ymax = Inf), alpha = 0.5) + 
  geom_line(data = snp, aes(date, pd)) +
  labs(title = "Figure 2: S&P 500 Price-to-Dividend Ratio",
       subtitle = "January 1973 - July 2018",
       caption = "Notes: The solid
line is the price-to-dividend ratio and the shaded areas are the periods where
the PSY statistic exceeds its 95% bootstrapped critical value.", 
       x = "Year", y = "Ratio") 

###

# NASDAQ

library(quantmod)

getSymbols('^IXIC',src='yahoo')
NASDAQ <- IXIC
class(NASDAQ)

head(NASDAQ)
tail(NASDAQ)
price <- NASDAQ[,4]
plot(price)