# Load necessary libraries
packages <- c(
    "vars",
    "ggplot2",
    "readr",
    "lubridate",
    "dplyr"
)

for (package in packages) {
    if (!require(package, character.only = TRUE)) {
        install.packages(package)
        library(package, character.only = TRUE)
    }
}

# Load the data
url <- 'https://www.fbc.keio.ac.jp/~tyabu/keiryo/oil2.csv'
df <- read.csv(url, header = TRUE)
df$date <- dmy(paste0("01-", df$date))
df <- df %>% arrange(date)
head(df)

# Fit the VAR model
var_model <- VAR(df[, -1], p = 4, type = "const")

# Determine the optimal lag order
lag_order <- VARselect(df[, -1], lag.max = 4, type = "const")
print(lag_order$selection)

# Estimate the VAR model
var_result <- VAR(df[, -1], p = lag_order$selection["AIC(n)"], type="const")
summary(var_result)

# Plot the impulse response functions
irf_oil <- irf(var_result, impulse = "oil", response = "price", n.ahead = 20, boot = TRUE, ci = 0.95)
plot(irf_oil)

irf_output <- irf(var_result, impulse = "output", response = "price", n.ahead = 20, boot = TRUE, ci = 0.95)
plot(irf_output)

irf_price <- irf(var_result, impulse = "price", response = "price", n.ahead = 20, boot = TRUE, ci = 0.95)
plot(irf_price)

# Extract IRF values
irf_values <- irf$irf

# Create a data frame for IRF values
irf_df <- data.frame(
    Period = 1:21,
    Oil = irf_oil['irf'],
    Oil_Lower = irf_oil['Lower'],
    Oil_Upper = irf_oil['Upper'],
    Output = irf_output['irf'],
    Output_Lower = irf_output['Lower'],
    Output_Upper = irf_output['Upper'],
    Price = irf_price['irf'],
    Price_Lower = irf_price['Lower'],
    Price_Upper = irf_price['Upper']
)

# Print the IRF table
# Rename the columns
colnames(irf_df) <- c(
    "Period",
    "Oil_IRF", "Oil_Lower_CI", "Oil_Upper_CI",
    "Output_IRF", "Output_Lower_CI", "Output_Upper_CI",
    "Price_IRF", "Price_Lower_CI", "Price_Upper_CI"
)

# Print the IRF table
print(irf_df)