Share


You said:

The second report should consist of an application of the time series econometrics learnt in Modules 4-5 to your data. Full submission should consist of a document containing the introduction, literature review, methods and data and preliminary results; the file containing your data and your computer codes. Please present the results (tables, plots, etc. in your report).

USL is United States 12 Month Oil Fund LP.
FUT is Future oil Prices.
My data:
> str(merged_data)
tibble [3,148 × 3] (S3: tbl_df/tbl/data.frame)
 $ Dates: POSIXct[1:3148], format: "2019-12-30" ...
 $ Fut  : num [1:3148] 61.7 61.7 61.7 61.1 61.1 ...
 $ USL  : num [1:3148] 23.2 23.2 23.2 23 23 ...
 - attr(*, "na.action")= 'omit' Named int [1:3] 3149 3150 3151
  ..- attr(*, "names")= chr [1:3] "3149" "3150" "3151"
> summary(merged_data)
     Dates                             Fut              USL       
 Min.   :2007-12-06 00:00:00.00   Min.   : 26.21   Min.   :13.60  
 1st Qu.:2010-12-12 06:00:00.00   1st Qu.: 52.77   1st Qu.:21.31  
 Median :2013-12-17 12:00:00.00   Median : 72.12   Median :36.84  
 Mean   :2013-12-17 07:11:49.01   Mean   : 73.98   Mean   :33.90  
 3rd Qu.:2016-12-22 06:00:00.00   3rd Qu.: 94.27   3rd Qu.:42.36  
 Max.   :2019-12-30 00:00:00.00   Max.   :145.29   Max.   :87.85  

My code:

# Load necessary libraries
library(readxl)  # For reading Excel files
library(dplyr)   # For data manipulation
library(tidyr)   # For tidying data
library(ggplot2) # For data visualization
library(MTS)     # For time series analysis
library(data.table) # For data handling
library(urca)    # For unit root tests
library(vars)    # For vector autoregression (VAR) modeling
library(car)     # For linear hypothesis testing
library(MCMCpack) # For Markov Chain Monte Carlo (MCMC) computations
library(maxLik)  # For maximum likelihood estimation
library(mvtnorm) # For multivariate normal distribution handling
library(rugarch) # For GARCH modeling

# Set working directory
setwd("C:/Users/abelma2/Desktop/R_Projects/ACE_527_Price_Discovery")

# Load datasets
Fut <- read_excel("OILNF_US0_NEW.xlsx", sheet = "Fut", col_types = c("date", "numeric"))
USL <- read_excel("OILNF_US0_NEW.xlsx", sheet = "USL", col_types = c("date", "numeric"))

# Merge datasets by date and remove rows with missing values
merged_data <- Fut %>%
  full_join(USL, by = "Dates") %>% 
  na.omit()

# Plot merged data
par(mfrow = c(2, 1), mar = c(3, 4, 1, 1))  # Adjust plotting parameters
names <- c("Fut", "USL")
for (name in names) {
  plot(merged_data$Dates, merged_data[[name]], type = "l", 
       ylab = "", 
       main = paste("Plot of", name))
}

# Transform data to logarithmic scale and add time column
log_data <- Fut %>%
  arrange(Dates) %>% 
  full_join(USL, by = "Dates") %>%
  na.omit() %>%
  mutate(
    log_Fut = log(Fut),
    log_USL = log(USL),
    Time = row_number()  # Add an incremental time column
  ) %>%
  select(Time, Dates, log_Fut, log_USL)

# Plot log-transformed data
par(mfrow = c(2, 1), mar = c(3, 4, 1, 1))  # Adjust plotting parameters
log_names <- c("log_Fut", "log_USL")
for (name in log_names) {
  plot(log_data$Dates, log_data[[name]], type = "l", 
       ylab = "", 
       main = paste("Plot of", name))
}

# Compute first differences of log-transformed data
diff_log_data <- log_data %>%
  arrange(Dates) %>%
  mutate(
    diff_log_Fut = c(NA, diff(log_Fut)),
    diff_log_USL = c(NA, diff(log_USL))
  ) %>%
  drop_na() %>%
  select(Time, Dates, diff_log_Fut, diff_log_USL)

# Plot first differences
par(mfrow = c(2, 1), mar = c(3, 4, 1, 1))  # Adjust plotting parameters
diff_log_names <- c("diff_log_Fut", "diff_log_USL")
for (name in diff_log_names) {
  plot(diff_log_data$Dates, diff_log_data[[name]], type = "l", 
       ylab = "", 
       main = paste("Plot of", name), 
       xlab = "Dates")
}

# Compute first differences
diff_log_dataframe <- diff_log_data %>%
  select(Dates, diff_log_Fut, diff_log_USL) %>% 
  as.data.frame()


# Conduct Dickey-Fuller tests on log-transformed data
dicky_fuller.trend <- lapply(log_data[, c("log_Fut", "log_USL")], 
                             function(x) ur.df(as.vector(x), type = "trend", selectlags = "AIC"))
dicky_fuller.constant <- lapply(log_data[, c("log_Fut", "log_USL")], 
                                function(x) ur.df(as.vector(x), type = "drift", selectlags = "AIC"))
dicky_fuller.none <- lapply(log_data[, c("log_Fut", "log_USL")], 
                            function(x) ur.df(as.vector(x), type = "none", selectlags = "AIC"))

# Print summaries of Dickey-Fuller tests
lapply(dicky_fuller.trend, summary)
lapply(dicky_fuller.constant, summary)
lapply(dicky_fuller.none, summary)

# Conduct Dickey-Fuller tests on first-difference data
diff_df.trend <- lapply(diff_log_data[, c("diff_log_Fut", "diff_log_USL")], 
                        function(x) ur.df(as.vector(x), type = "trend", selectlags = "AIC"))
diff_df.constant <- lapply(diff_log_data[, c("diff_log_Fut", "diff_log_USL")], 
                           function(x) ur.df(as.vector(x), type = "drift", selectlags = "AIC"))
diff_df.none <- lapply(diff_log_data[, c("diff_log_Fut", "diff_log_USL")], 
                       function(x) ur.df(as.vector(x), type = "none", selectlags = "AIC"))

# Print summaries of Dickey-Fuller tests on differences
lapply(diff_df.trend, summary)
lapply(diff_df.constant, summary)
lapply(diff_df.none, summary)

# Identify optimal lag order for VAR model
optimal <- VARorder(diff_log_dataframe[, c("diff_log_Fut", "diff_log_USL")], maxp = 100)

detach("package:MTS", unload = TRUE)

# Fit VAR model
var.m1 <- VAR(diff_log_dataframe[, c("diff_log_Fut", "diff_log_USL")], p = 6, type = "const")
summary(var.m1)

var.m2 <- VAR(diff_log_dataframe[, c("diff_log_USL", "diff_log_Fut")], p = 6, type = "const")
summary(var.m2)

# Conduct hypothesis tests on VAR model coefficients
vars <- c("diff_log_Fut", "diff_log_USL")
names <- c("diff_log_Fut.l", "diff_log_USL.l")
for (n in names) {
  for (var in vars) {
    hypothesis_matrix <- paste0(n, 1:6)  # Create hypothesis matrix
    result <- linearHypothesis(
      model = var.m1$varresult[[var]], 
      hypothesis.matrix = hypothesis_matrix, 
      rhs = rep(0, 6)
    )
    print(result)
  }
}

# Conduct serial correlation test on VAR model
serial.test(var.m1, lags.pt = 20, type = "PT.asymptotic")
roots(var.m1)

# Compute AIC and SBC for VAR model
res <- residuals(var.m1)
N <- var.m1$obs
nreg <- var.m1$K * (var.m1$p * var.m1$K + 1)
det.cov.m1 <- det(crossprod(res) / N)
aic.var.m1 <- N * log(det.cov.m1) + 2 * nreg
sbc.var.m1 <- N * log(det.cov.m1) + log(N) * nreg
print(paste("AIC.var.m1 =", sprintf("%.3f", aic.var.m1), 
            "SBC.var.m1 =", sprintf("%.3f", sbc.var.m1)))

# Compute and plot Forecast Error Variance Decomposition (FEVD) for the VAR models
fevd.m1 <- fevd(var.m1, n.ahead = 30)
plot(fevd.m1, 
     main = "FEVD for VAR Model 1", 
     col = c("blue", "red"),  # Color for the two series
     lwd = 2,                # Line width
     xlab = "Forecast Horizon", 
     ylab = "Variance Decomposition", 
     col.axis = "black",     # Color for axis labels
     col.lab = "black",      # Color for axis titles
     col.main = "black"      # Color for the main title
)

fevd.m2 <- fevd(var.m2, n.ahead = 30)
plot(fevd.m2, 
     main = "FEVD for VAR Model 2", 
     col = c("red", "blue"),  # Color for the two series
     lwd = 2,                    # Line width
     xlab = "Forecast Horizon", 
     ylab = "Variance Decomposition", 
     col.axis = "black",         # Color for axis labels
     col.lab = "black",          # Color for axis titles
     col.main = "black"          # Color for the main title
)

# Compute Impulse Response Function (IRF) for VAR model 1
irf.m1 <- irf(var.m1, impulse = NULL, response = NULL, n.ahead = 30, 
              ortho = TRUE, boot = TRUE, runs = 100)

# Access and plot IRF results for the two series
irf.m1$irf$diff_log_Fut
irf.m1$irf$diff_log_USL

# Prepare and plot the impulse response data
plot.data <- rbindlist(
  lapply(
    irf.m1[c("irf", "Lower", "Upper")],
    function(dt) {
      temp <- rbindlist(
        lapply(dt, function(x) as.data.table(x)[, Step := 0:(.N-1)]), idcol = "Shock"
      )
      temp <- melt(temp, id.vars = c("Shock", "Step"), variable.name = "Response")
    }
  ),
  idcol = "data.type"
)

# Reorder 'Shock' and 'Response' factors in the desired order
plot.data[, Shock := factor(Shock, levels = c("diff_log_Fut", "diff_log_USL"))]
plot.data[, Response := factor(Response, levels = c("diff_log_Fut", "diff_log_USL"))]

# Create the plot with reordered factors
ggplot(data = plot.data,
       aes(x = Step, y = value, linetype = data.type)) + 
  geom_line() + 
  facet_grid(rows = vars(Response), cols = vars(Shock), scales = "free_y") +
  geom_hline(yintercept = 0, colour = "red") + 
  theme_bw() + 
  theme(legend.position = "none") +
  ylab("") + 
  ggtitle("Responses (in rows) to Shocks (in columns)")

# Plot all residuals on the same graph
matplot(res, type = "l", col = 1:5, lty = 1, xlab = "Index", ylab = "Residuals", main = "Residuals of VAR Model")

# Add a legend to the residuals plot
legend("topright", legend = colnames(res), col = 1:5, lty = 1)

# Fit GARCH(1,1) models to residuals for volatility modeling
garch11 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                      mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
                      distribution.model = "norm")

garch.fit.FUT <- ugarchfit(res[, 1], spec = garch11)
garch.fit.USL <- ugarchfit(res[, 2], spec = garch11)

# Add GARCH conditional variances to the dataset for further analysis
aligned_diff_log_dataframe <- diff_log_dataframe[-c(1:(nrow(diff_log_dataframe) - length(coredata(sigma(garch.fit.FUT))^2))), ]
aligned_diff_log_dataframe$f_t <- coredata(sigma(garch.fit.FUT))^2
aligned_diff_log_dataframe$U_t <- coredata(sigma(garch.fit.USL))^2

# Plot GARCH(1,1) conditional variances for both series
par(mfrow = c(2, 1))
plot(aligned_diff_log_dataframe$Dates, aligned_diff_log_dataframe$f_t, type = "l",
     main = "GARCH(1,1) conditional variance of FUT", ylab = "")
plot(aligned_diff_log_dataframe$Dates, aligned_diff_log_dataframe$U_t, type = "l",
     main = "GARCH(1,1) conditional variance of USL", ylab = "")

# Plot standardized residuals from the GARCH models
stdres_Fut <- coredata(residuals(garch.fit.FUT, standardize = TRUE))
stdres_USL <- coredata(residuals(garch.fit.USL, standardize = TRUE))
plot(stdres_Fut, main = "Standardized Residuals (FUT)")
plot(stdres_USL, main = "Standardized Residuals (USL)")