suppressMessages(suppressWarnings({
  library(robustbase)
  library("Hmisc")
  library("plm")
  library("lmtest")
  library("stargazer")
  library("ivreg")
  library("systemfit")
  library("Hmisc")
  library(ggplot2)
  library("pscl")
  library(dplyr)
  library("stringr")
  library(lmtest)
  library(censReg)
  library(mfx)
  library(corrplot)
  library(pscl)
  library(ROCR)
  library(caret)
  library(pROC)
  library(corrplot)
  library(dplyr)
  library(tidyverse)
  library("mlogit")
  library(patchwork) # To display 2 charts together
  library(scales)
  library("reshape2")
  library(ggplot2)
  library(tseries)
  library(seasonal)
}))

#Loading the data, converting them to time series, calculating logarithm and performing adf test for GDP

#loading the data
gdp_q_data <- read.csv("GDP.csv", header=TRUE, sep=",", dec=".")
head(gdp_q_data)
# data frame
gdp_q_data <- data.frame(
  QYEAR = gdp_q_data$Year,
  GDP = gdp_q_data$GDP)
# Parse the YEAR column to extract year and quarter
gdp_q_data$Year <- as.numeric(sub(".*Q(\\d{4})$", "\\1", gdp_q_data$QYEAR))
gdp_q_data$Quarter <- as.numeric(sub("Q.*", "", gdp_q_data$QYEAR))
# Print the modified data frame to verify parsing
print(gdp_q_data)

# Create a time series object for the GDP column
# Start year: min(gdp_q_data$Year)
# Frequency: 4 (quarterly data)
gdp_ts <- ts(gdp_q_data$GDP, start = c(min(gdp_q_data$Year), min(gdp_q_data$Quarter)), frequency = 4)
# Print the time series object
print(gdp_ts)
# Plot the time series
plot(gdp_ts, main = "Quarterly GDP Data", ylab = "GDP", xlab = "Year")

# Calculate logarithm of GDP and add it as a new variable
log_gdp_ts <- log(gdp_ts)
# Print the updated data
print(log_gdp_ts)
# Plot log GDP
plot(log_gdp_ts, main = "Quarterly Log GDP Data", ylab = "Log GDP", xlab = "Year")
library(tseries)
# Perform ADF test on log GDP
adf_result_log_gdp_ts <- adf.test(log_gdp_ts)
# Print the ADF test result
print(adf_result_log_gdp_ts) #not stationary



#Loading the data, converting them to time series, calculating logarithm and performing adf test for M1

#loading the data for M1
M1_q_data <- read.csv("M1q.csv", header=TRUE, sep=",", dec=".")
head(M1_q_data)
# data frame
M1_q_data <- data.frame(
  QYEAR = M1_q_data$Year,
  M1 = M1_q_data$M1)
# Parse the YEAR column to extract year and quarter
M1_q_data$Year <- as.numeric(sub(".*Q(\\d{4})$", "\\1", M1_q_data$QYEAR))
M1_q_data$Quarter <- as.numeric(sub("Q.*", "", M1_q_data$QYEAR))
# Print the modified data frame to verify parsing
print(M1_q_data)

# Create a time series object for the M1 column
# Start year: min(M1_q_data$Year)
# Frequency: 4 (quarterly data)
m1_ts <- ts(M1_q_data$M1, start = c(min(M1_q_data$Year), min(M1_q_data$Quarter)), frequency = 4)
# Print the time series object
print(m1_ts)
# Plot the time series
plot(m1_ts, main = "Quarterly M1 Data", ylab = "M1", xlab = "Year")

# Calculate logarithm of M1 and add it as a new variable
log_m1_ts <- log(m1_ts)
# Print the updated data
print(log_m1_ts)
# Plot log M1
plot(log_m1_ts, main = "Quarterly Log M1 Data", ylab = "Log M1", xlab = "Year")

# Perform ADF test on log M1
adf_result_log_m1_ts <- adf.test(log_m1_ts)
# Print the ADF test result
print(adf_result_log_m1_ts) #not stationary



#Loading the data, converting them to time series, calculating logarithm and performing adf test for M2

#loading the data for M2
M2_q_data <- read.csv("M2q.csv", header=TRUE, sep=",", dec=".")
head(M2_q_data)
# data frame
M2_q_data <- data.frame(
  QYEAR = M2_q_data$Year,
  M2 = M2_q_data$M2)
# Parse the YEAR column to extract year and quarter
M2_q_data$Year <- as.numeric(sub(".*Q(\\d{4})$", "\\1", M2_q_data$QYEAR))
M2_q_data$Quarter <- as.numeric(sub("Q.*", "", M2_q_data$QYEAR))
# Print the modified data frame to verify parsing
print(M2_q_data)

# Create a time series object for the M2 column
# Start year: min(M2_q_data$Year)
# Frequency: 4 (quarterly data)
m2_ts <- ts(M2_q_data$M2, start = c(min(M2_q_data$Year), min(M2_q_data$Quarter)), frequency = 4)
# Print the time series object
print(m2_ts)
# Plot the time series
plot(m2_ts, main = "Quarterly M2 Data", ylab = "M2", xlab = "Year")

# Calculate logarithm of M2 and add it as a new variable
log_m2_ts <- log(m2_ts)
# Print the updated data
print(log_m2_ts)
# Plot log M1
plot(log_m2_ts, main = "Quarterly Log M2 Data", ylab = "Log M2", xlab = "Year")

# Perform ADF test on log M2
adf_result_log_m2_ts <- adf.test(log_m2_ts)
# Print the ADF test result
print(adf_result_log_m2_ts) #not stationary



#Loading the data, converting them to time series, calculating logarithm and performing adf test for M3

#loading the data for M3
M3_q_data <- read.csv("M3q.csv", header=TRUE, sep=",", dec=".")
head(M3_q_data)
#data frame
M3_q_data <- data.frame(
  QYEAR = M3_q_data$Year,
  M3 = M3_q_data$M3)
# Parse the YEAR column to extract year and quarter
M3_q_data$Year <- as.numeric(sub(".*Q(\\d{4})$", "\\1", M3_q_data$QYEAR))
M3_q_data$Quarter <- as.numeric(sub("Q.*", "", M3_q_data$QYEAR))
# Print the modified data frame to verify parsing
print(M3_q_data)

# Create a time series object for the M3 column
# Start year: min(M3_q_data$Year)
# Frequency: 4 (quarterly data)
m3_ts <- ts(M3_q_data$M3, start = c(min(M3_q_data$Year), min(M3_q_data$Quarter)), frequency = 4)

# Print the time series object
print(m3_ts)
# Plot the time series
plot(m3_ts, main = "Quarterly M3 Data", ylab = "M3", xlab = "Year")

# Calculate logarithm of M3 and add it as a new variable
log_m3_ts <- log(m3_ts)
# Print the updated data
print(log_m3_ts)
# Plot log M3
plot(log_m3_ts, main = "Quarterly Log M3 Data", ylab = "Log M3", xlab = "Year")

# Perform ADF test on log M3
adf_result_log_m3_ts <- adf.test(log_m3_ts)
# Print the ADF test result
print(adf_result_log_m3_ts) #not stationary



# 1st diff log
d1_log_GDP <- diff(log_gdp_ts)
plot.ts(d1_log_GDP, xlab= 'Time', ylab='d1 Log GDP')

d1_log_M1<- diff(log_m1_ts)
plot.ts(d1_log_M1, xlab= 'Time', ylab='d1 Log M1')

d1_log_M2 <- diff(log_m2_ts)
plot.ts(d1_log_M2,xlab= 'Time', ylab='d1 Log M2')

d1_log_M3 <- diff(log_m3_ts)
plot.ts(d1_log_M3,xlab= 'Time', ylab='d1 Log M3')



# ADF test 1st diff log
adf.test(d1_log_GDP) #I(1) (stationary)
adf.test(d1_log_M1) #at least I(1) (not stationary)
adf.test(d1_log_M2) #I(1) (stationary)
adf.test(d1_log_M3) #I(1) (stationary)

# 2nd diff
d2_log_M1 <- diff(d1_log_M1)
plot.ts(d2_log_M1)
adf.test(d2_log_M1) #I(2)


#############################################################################

# Define a vector of lag values for GDP~M1
lag_values <- 1:80

# Initialize empty vectors to store coefficient estimates and confidence intervals
coefficients <- numeric(length(lag_values))
lower_ci <- numeric(length(lag_values))
upper_ci <- numeric(length(lag_values))

# Loop through each lag value
for (i in 1:length(lag_values)) {
  # Calculate the differences with lag for log GDP and log M1 money supply
  y <- diff(log_gdp_ts, lag = lag_values[i])  # Difference in log GDP with lag
  x <- diff(log_m1_ts, lag = lag_values[i])   # Difference in log M1 money supply with lag
 
  # Estimate the linear regression model
  ols <- lm(y ~ x)
 
  # Store coefficient estimate for the lag
  coefficients[i] <- coef(ols)[2]  # Assuming x is the second variable in the model
 
  # Compute confidence intervals for the coefficient estimate
  ci <- confint(ols)[2, ]  # Confidence interval for the coefficient of x
  lower_ci[i] <- ci[1]
  upper_ci[i] <- ci[2]
}

# Plot coefficient estimates with confidence intervals against lag values
plot(lag_values, coefficients, type = "b", pch = 19,
     xlab = "Lag (j)", ylab = "Coefficient Estimate",
     main = "M1")
lines(lag_values, upper_ci, col = "red", lty = 2)  # Upper confidence interval
lines(lag_values, lower_ci, col = "blue", lty = 2)  # Lower confidence interval
lines(rep(0,length(lag_values)), col="green")
legend("topleft", legend = c("Coefficient Estimate", "Upper CI", "Lower CI"),
       col = c("black", "red", "blue"), lty = c(1, 2, 2), pch = c(19, NA, NA))



# Define a data frame to store lag values, coefficient estimates, and confidence intervals
results <- data.frame(
  lag = lag_values,
  coefficient = coefficients,
  lower_ci = lower_ci,
  upper_ci = upper_ci
)

ggplot(results, aes(x = lag, y = coefficient)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2, color = "blue") +
  geom_line(color = "black", linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "solid", color = "red") +  # Add horizontal line at y = 0
  labs(x = "Lag (j)", y = "Coefficient Estimate",
       title = "Coefficient Estimates with Confidence Intervals") +
  theme_minimal()

#############################

# Define a vector of lag values for GDP~M2
lag_values <- 1:80

# Initialize empty vectors to store coefficient estimates and confidence intervals
coefficients <- numeric(length(lag_values))
lower_ci <- numeric(length(lag_values))
upper_ci <- numeric(length(lag_values))

# Loop through each lag value
for (i in 1:length(lag_values)) {
  # Calculate the differences with lag for log GDP and log M2 money supply
  y <- diff(log_gdp_ts, lag = lag_values[i])  # Difference in log GDP with lag
  x <- diff(log_m2_ts, lag = lag_values[i])   # Difference in log M2 money supply with lag
 
  # Estimate the linear regression model
  ols <- lm(y ~ x)
 
  # Store coefficient estimate for the lag
  coefficients[i] <- coef(ols)[2]  # Assuming x is the second variable in the model
 
  # Compute confidence intervals for the coefficient estimate
  ci <- confint(ols)[2, ]  # Confidence interval for the coefficient of x
  lower_ci[i] <- ci[1]
  upper_ci[i] <- ci[2]
}

# Plot coefficient estimates with confidence intervals against lag values
plot(lag_values, coefficients, type = "b", pch = 19,
     xlab = "Lag (j)", ylab = "Coefficient Estimate",
     main = "M2")
lines(lag_values, upper_ci, col = "red", lty = 2)  # Upper confidence interval
lines(lag_values, lower_ci, col = "blue", lty = 2)  # Lower confidence interval
lines(rep(0,length(lag_values)), col="green")
legend("topleft", legend = c("Coefficient Estimate", "Upper CI", "Lower CI"),
       col = c("black", "red", "blue"), lty = c(1, 2, 2), pch = c(19, NA, NA))



# Define a data frame to store lag values, coefficient estimates, and confidence intervals
results <- data.frame(
  lag = lag_values,
  coefficient = coefficients,
  lower_ci = lower_ci,
  upper_ci = upper_ci
)

ggplot(results, aes(x = lag, y = coefficient)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2, color = "blue") +
  geom_line(color = "black", linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "solid", color = "red") +  # Add horizontal line at y = 0
  labs(x = "Lag (j)", y = "Coefficient Estimate",
       title = "Coefficient Estimates with Confidence Intervals") +
  theme_minimal()

#############################

# Define a vector of lag values for GDP~M3
lag_values <- 1:80

# Initialize empty vectors to store coefficient estimates and confidence intervals
coefficients <- numeric(length(lag_values))
lower_ci <- numeric(length(lag_values))
upper_ci <- numeric(length(lag_values))

# Loop through each lag value
for (i in 1:length(lag_values)) {
  # Calculate the differences with lag for log GDP and log M1 money supply
  y <- diff(log_gdp_ts, lag = lag_values[i])  # Difference in log GDP with lag
  x <- diff(log_m3_ts, lag = lag_values[i])   # Difference in log M3 money supply with lag
 
  # Estimate the linear regression model
  ols <- lm(y ~ x)
 
  # Store coefficient estimate for the lag
  coefficients[i] <- coef(ols)[2]  # Assuming x is the second variable in the model
 
  # Compute confidence intervals for the coefficient estimate
  ci <- confint(ols)[2, ]  # Confidence interval for the coefficient of x
  lower_ci[i] <- ci[1]
  upper_ci[i] <- ci[2]
}

# Plot coefficient estimates with confidence intervals against lag values
plot(lag_values, coefficients, type = "b", pch = 19,
     xlab = "Lag (j)", ylab = "Coefficient Estimate",
     main = "M3")
lines(lag_values, upper_ci, col = "red", lty = 2)  # Upper confidence interval
lines(lag_values, lower_ci, col = "blue", lty = 2)  # Lower confidence interval
lines(rep(0,length(lag_values)), col="green")
legend("topleft", legend = c("Coefficient Estimate", "Upper CI", "Lower CI"),
       col = c("black", "red", "blue"), lty = c(1, 2, 2), pch = c(19, NA, NA))



# Define a data frame to store lag values, coefficient estimates, and confidence intervals
results <- data.frame(
  lag = lag_values,
  coefficient = coefficients,
  lower_ci = lower_ci,
  upper_ci = upper_ci
)

ggplot(results, aes(x = lag, y = coefficient)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2, color = "blue") +
  geom_line(color = "black", linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "solid", color = "red") +  # Add horizontal line at y = 0
  labs(x = "Lag (j)", y = "Coefficient Estimate",
       title = "Coefficient Estimates with Confidence Intervals") +
  theme_minimal()



# Define a vector of lag values for GDP~d_M1 (superneutrality)
lag_values <- 1:80

# Initialize empty vectors to store coefficient estimates and confidence intervals
coefficients <- numeric(length(lag_values))
lower_ci <- numeric(length(lag_values))
upper_ci <- numeric(length(lag_values))

# Loop through each lag value
for (i in 1:length(lag_values)) {
  # Calculate the differences with lag for log GDP and diff log M1 money supply
  y <- diff(log_gdp_ts, lag = lag_values[i])  # Difference in log GDP with lag
  x <- diff(d1_log_M1, lag = lag_values[i])   # Difference in diff log M1 money supply with lag
 
  # Ensure x and y are of the same length (already first difference M1)
  min_length <- min(length(y), length(x))
  y <- y[1:min_length]
  x <- x[1:min_length]
 
  # Estimate the linear regression model
  ols <- lm(y ~ x)
 
  # Store coefficient estimate for the lag
  coefficients[i] <- coef(ols)[2]  # Assuming x is the second variable in the model
 
  # Compute confidence intervals for the coefficient estimate
  ci <- confint(ols)[2, ]  # Confidence interval for the coefficient of x
  lower_ci[i] <- ci[1]
  upper_ci[i] <- ci[2]
}

# Plot coefficient estimates with confidence intervals against lag values
plot(lag_values, coefficients, type = "b", pch = 19,
     xlab = "Lag (j)", ylab = "Coefficient Estimate",
     main = "M1")
lines(lag_values, upper_ci, col = "red", lty = 2)  # Upper confidence interval
lines(lag_values, lower_ci, col = "blue", lty = 2)  # Lower confidence interval
lines(rep(0,length(lag_values)), col="green")
legend("topleft", legend = c("Coefficient Estimate", "Upper CI", "Lower CI"),
       col = c("black", "red", "blue"), lty = c(1, 2, 2), pch = c(19, NA, NA))



# Define a data frame to store lag values, coefficient estimates, and confidence intervals
results <- data.frame(
  lag = lag_values,
  coefficient = coefficients,
  lower_ci = lower_ci,
  upper_ci = upper_ci
)

ggplot(results, aes(x = lag, y = coefficient)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2, color = "blue") +
  geom_line(color = "black", linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "solid", color = "red") +  # Add horizontal line at y = 0
  labs(x = "Lag (j)", y = "Coefficient Estimate",
       title = "Coefficient Estimates with Confidence Intervals") +
  theme_minimal()