if (!requireNamespace("readxl", quietly = TRUE)) {
  install.packages("readxl")
}

if (!requireNamespace("rmarkdown", quietly = TRUE)) {
  install.packages("rmarkdown")
}

# Load the readxl package
library(readxl)

# Specify the file path
file_path <- "Final.xlsx"

# Read the Excel file into a data frame
data <- read_excel(file_path)
## New names:
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
data <- subset(data, select = 1:(ncol(data) - 4))

data <- data[, !(names(data) %in% c("gas_f", "gas_s"))]

# Assuming 'data' is your data frame
data$s <- log(data$soy_s)
data$f <- log(data$soy_f)

head(data)
## # A tibble: 6 × 5
##   Date                soy_f soy_s     s     f
##   <dttm>              <dbl> <dbl> <dbl> <dbl>
## 1 2006-06-13 00:00:00  586.   562  6.33  6.37
## 2 2006-06-20 00:00:00  600.   579  6.36  6.40
## 3 2006-06-27 00:00:00  580.   560  6.33  6.36
## 4 2006-07-04 00:00:00  599.   574  6.35  6.40
## 5 2006-07-11 00:00:00  606.   577  6.36  6.41
## 6 2006-07-18 00:00:00  603    577  6.36  6.40
if (!requireNamespace("ggplot2", quietly = TRUE)) {
  install.packages("ggplot2")
}

# Assuming 'data' is your data frame and you have the ggplot2 library installed
library(ggplot2)
library(scales)  # Required for date_format function

# Convert 'Date' column to Date class
data$Date <- as.Date(data$Date)

# Create a ggplot with custom aesthetics
my_plot <- ggplot(data, aes(x = Date)) +
  geom_line(aes(y = s, color = "Spot"), size = 0.5) +
  geom_line(aes(y = f, color = "Forward"), size = 1, linetype = "dashed") +
  labs(x = "Year", y = "Natural Logarithm of Price") +
  scale_color_manual(values = c("Spot" = "steelblue", "Forward" = "coral")) +
  
  # Scale x-axis with breaks for years and minor breaks for months
  scale_x_date(
    date_breaks = "1 year",  # Set major breaks for years
    date_labels = "%Y",      # Format major break labels as years
    labels = date_format("%b", tz = "UTC"),  # Format all labels as months
    expand = c(0, 0)
  ) +
  
  # Add a theme for better appearance
  theme_minimal() +
  theme(
    legend.position = "bottom",  # Move the legend to the bottom
    panel.grid.major = element_blank(),  # Remove major gridlines
    panel.grid.minor = element_blank(),  # Remove minor gridlines
    axis.line = element_line(color = "black"),  # Set color of axis lines
    legend.box.background = element_rect(color = "black"),  # Set legend box color
    legend.text = element_text(size = 10),  # Adjust legend text size
    legend.title = element_blank(),  # Remove legend title
    axis.text = element_text(size = 10),  # Adjust axis text size
    axis.title = element_text(size = 10, face = "bold")  # Adjust axis title size and style
  )
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
my_plot

# Save the plot as a PNG file
ggsave("C:/Users/ehsan/OneDrive - University of Essex/Courses/Econometrics/Coursework2/R/Question1.png", plot = my_plot, width = 10, height = 6, dpi = 300)

if (!requireNamespace("tseries", quietly = TRUE)) {
  install.packages("tseries")
}
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)

# Install and load the urca package if not already installed
if (!requireNamespace("urca", quietly = TRUE)) {
  install.packages("urca")
}
library(urca)

# Assuming 'data' is your data frame
data$Date <- as.Date(data$Date)

# Extract 's' and 'f' columns as time series
s_time_series <- ts(data$s, frequency = 52)  # Assuming weekly data
f_time_series <- ts(data$f, frequency = 52)  # Assuming weekly data

# Perform Augmented Dickey-Fuller test with automatic lag selection for 's'
adf_test_s <- ur.df(s_time_series, type = "drift", selectlags = "BIC")

# Perform Augmented Dickey-Fuller test with automatic lag selection for 'f'
adf_test_f <- ur.df(f_time_series, type = "drift", selectlags = "BIC")

# Print the results
cat("ADF Test Results for 's':\n")
## ADF Test Results for 's':
summary(adf_test_s)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.158558 -0.018027  0.001761  0.023835  0.129068 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.113027   0.041314   2.736  0.00640 **
## z.lag.1     -0.016099   0.005924  -2.717  0.00676 **
## z.diff.lag   0.067126   0.040115   1.673  0.09477 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03647 on 612 degrees of freedom
## Multiple R-squared:  0.01584,    Adjusted R-squared:  0.01263 
## F-statistic: 4.926 on 2 and 612 DF,  p-value: 0.007545
## 
## 
## Value of test-statistic is: -2.7175 3.8507 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
cat("\nADF Test Results for 'f':\n")
## 
## ADF Test Results for 'f':
summary(adf_test_f)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.159983 -0.018101  0.001826  0.022947  0.102281 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.115848   0.042074   2.753  0.00607 **
## z.lag.1     -0.016465   0.006022  -2.734  0.00644 **
## z.diff.lag   0.026860   0.040188   0.668  0.50415   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0349 on 612 degrees of freedom
## Multiple R-squared:  0.01257,    Adjusted R-squared:  0.009347 
## F-statistic: 3.897 on 2 and 612 DF,  p-value: 0.02082
## 
## 
## Value of test-statistic is: -2.734 3.933 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
# Take the first differences of 's' and 'f'
s_diff <- diff(s_time_series)
f_diff <- diff(f_time_series)

# Perform Augmented Dickey-Fuller test for first differences of 's' with only an intercept
adf_test_s_diff <- ur.df(s_diff, type = "none", selectlags = "BIC")

# Perform Augmented Dickey-Fuller test for first differences of 'f' with only an intercept
adf_test_f_diff <- ur.df(f_diff, type = "none", selectlags = "BIC")

# Print the results
cat("ADF Test Results for first differences of 's':\n")
## ADF Test Results for first differences of 's':
summary(adf_test_s_diff)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.153844 -0.018249  0.003462  0.022235  0.134815 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.88027    0.05517 -15.955   <2e-16 ***
## z.diff.lag -0.05869    0.04032  -1.456    0.146    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03661 on 612 degrees of freedom
## Multiple R-squared:  0.4697, Adjusted R-squared:  0.468 
## F-statistic:   271 on 2 and 612 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.9548 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62
cat("\nADF Test Results for first differences of 'f':\n")
## 
## ADF Test Results for first differences of 'f':
summary(adf_test_f_diff)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15778 -0.01747  0.00348  0.02297  0.09580 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.92616    0.05637 -16.429   <2e-16 ***
## z.diff.lag -0.05063    0.04034  -1.255     0.21    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03505 on 612 degrees of freedom
## Multiple R-squared:  0.4894, Adjusted R-squared:  0.4877 
## F-statistic: 293.3 on 2 and 612 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -16.4289 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62
# Cointegration ###################
# Step 1: Run a linear regression of 's' on 'f'
regression_step1 <- lm(s_time_series ~ f_time_series)

# Step 2: Obtain the residuals from the regression
residuals_step2 <- residuals(regression_step1)

# Step 3: Test for stationarity of residuals (Augmented Dickey-Fuller test)
adf_test_residuals <- ur.df(residuals_step2, type = "none", selectlags = "BIC")

# Print the results
cat("ADF Test Results for residuals:\n")
## ADF Test Results for residuals:
summary(adf_test_residuals)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.087745 -0.003462 -0.000429  0.003045  0.120310 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.10810    0.01751  -6.174 1.22e-09 ***
## z.diff.lag  0.13721    0.04004   3.427 0.000651 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01293 on 613 degrees of freedom
## Multiple R-squared:  0.06498,    Adjusted R-squared:  0.06193 
## F-statistic:  21.3 on 2 and 613 DF,  p-value: 1.138e-09
## 
## 
## Value of test-statistic is: -6.1736 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62
# ECM #############

# Step 1: Run a cointegrating regression: s ~ f
cointegration_regression <- lm(s_time_series ~ f_time_series)

# Step 2: Obtain the residuals from the cointegrating regression
residuals_cointegration <- residuals(cointegration_regression)

# Step 3: Calculate first differences of 's' and 'f'
s_diff <- diff(s_time_series)
f_diff <- diff(f_time_series)

# Step 4: Create lagged residuals (lagged u_hat_t)
lagged_residuals <- c(NA, residuals_cointegration)

# Trim the vectors to have the same length
#s_diff <- s_diff[-1]
#f_diff <- f_diff[-1]
lagged_residuals <- c(NA, head(residuals_cointegration, -2))

# Combine variables into a data frame
ecm_data <- data.frame(s_diff = s_diff, f_diff = f_diff, lagged_residuals = lagged_residuals)

# Step 5: Formulate and estimate the Error Correction Model (ECM)
ecm_model <- lm(s_diff ~ f_diff + lagged_residuals, data = ecm_data)

# Print the results
summary(ecm_model)
## 
## Call:
## lm(formula = s_diff ~ f_diff + lagged_residuals, data = ecm_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.092973 -0.003463 -0.000309  0.002871  0.112029 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.256e-05  5.083e-04   0.044    0.965    
## f_diff            9.701e-01  1.455e-02  66.689  < 2e-16 ***
## lagged_residuals -1.182e-01  1.671e-02  -7.074 4.13e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0126 on 612 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.8825, Adjusted R-squared:  0.8821 
## F-statistic:  2298 on 2 and 612 DF,  p-value: < 2.2e-16
# Create R ###########
# Install and load the dplyr package if not already installed
if (!requireNamespace("dplyr", quietly = TRUE)) {
  install.packages("dplyr")
}
library(dplyr)
## 
## Attaching package: 'dplyr'
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Calculate log returns
log_returns <- log(data$soy_s/lag(data$soy_s, 1))
head(log_returns)
## [1]          NA  0.02980063 -0.03336569  0.02469261  0.00521287  0.00000000
# Create a new data frame with log returns and dates
returns_data <- data.frame(Date = data$Date, Log_Returns = log_returns)

# Specify the width for the PNG file
#png("C:/Users/ehsan/OneDrive - University of Essex/Courses/Econometrics/Coursework2/R/returns.png", width = 800)  # Adjust the width as needed

# Plot the return series
plot(returns_data$Date, returns_data$Log_Returns, type = "l", col = "blue", xlab = "Date", ylab = "Returns", main = "Log Return Series")

# Save and close the PNG device
#dev.off()

# AR(1) #####

# Estimate AR(1) model
ar1_model <- arima(log_returns, order = c(1, 0, 0))

# Load the required packages
if (!requireNamespace("stargazer", quietly = TRUE)) {
  install.packages("stargazer")
}
library(stargazer)
## 
## Please cite as: 
## 
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
# Print the table
stargazer(ar1_model, type = "text")
## 
## =============================================
##                       Dependent variable:    
##                   ---------------------------
##                           log_returns        
## ---------------------------------------------
## ar1                          0.063           
##                             (0.040)          
##                                              
## intercept                    0.001           
##                             (0.002)          
##                                              
## ---------------------------------------------
## Observations                  616            
## Log Likelihood             1,163.661         
## sigma2                       0.001           
## Akaike Inf. Crit.         -2,321.322         
## =============================================
## Note:             *p<0.1; **p<0.05; ***p<0.01
# Extract residuals
residuals_ar1 <- residuals(ar1_model)

# Calculate squared residuals
squared_residuals <- residuals_ar1^2

# Specify the number of lags for the ARCH test
num_lags_arch_test <- 5

if (!requireNamespace("FinTS", quietly = TRUE)) {
  install.packages("FinTS")
}

library(FinTS)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Perform the ARCH test
arch_test_result <- ArchTest(squared_residuals, lags = num_lags_arch_test , demean = TRUE)

# Print the ARCH test results
print(arch_test_result)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  squared_residuals
## Chi-squared = 25.02, df = 5, p-value = 0.0001381
# GARCH #####

# Load the required packages
if (!requireNamespace("rugarch", quietly = TRUE)) {
  install.packages("rugarch")
}
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## 
## The following object is masked from 'package:stats':
## 
##     sigma
if (any(is.na(log_returns))) {
  log_returns <- na.omit(log_returns)
}

# Specify the model
spec <- ugarchspec(mean.model = list(armaOrder = c(1, 0)), 
                   variance.model = list(model = "sGARCH", garchOrder = c(1, 1)))

# Estimate the model
garch_model <- ugarchfit(spec, data = log_returns)

# Print the model summary
show(garch_model)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001339    0.001240  1.08008 0.280105
## ar1     0.041768    0.045314  0.92176 0.356652
## omega   0.000133    0.000041  3.28006 0.001038
## alpha1  0.237694    0.054141  4.39029 0.000011
## beta1   0.672790    0.061746 10.89600 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001339    0.001276  1.04945 0.293970
## ar1     0.041768    0.045078  0.92659 0.354141
## omega   0.000133    0.000064  2.07080 0.038377
## alpha1  0.237694    0.059697  3.98168 0.000068
## beta1   0.672790    0.080190  8.38990 0.000000
## 
## LogLikelihood : 1219.842 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.9443
## Bayes        -3.9084
## Shibata      -3.9444
## Hannan-Quinn -3.9303
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.00802  0.9286
## Lag[2*(p+q)+(p+q)-1][2]   0.15316  0.9996
## Lag[4*(p+q)+(p+q)-1][5]   1.82227  0.7618
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.002485  0.9602
## Lag[2*(p+q)+(p+q)-1][5]  0.852251  0.8921
## Lag[4*(p+q)+(p+q)-1][9]  2.567736  0.8276
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1424 0.500 2.000  0.7059
## ARCH Lag[5]    1.0734 1.440 1.667  0.7113
## ARCH Lag[7]    1.4753 2.315 1.543  0.8263
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.5321
## Individual Statistics:              
## mu     0.29200
## ar1    0.02046
## omega  0.35968
## alpha1 0.36693
## beta1  0.58605
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.28 1.47 1.88
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.3175 0.7510    
## Negative Sign Bias  0.4570 0.6479    
## Positive Sign Bias  0.1139 0.9094    
## Joint Effect        0.2240 0.9736    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     31.08      0.03958
## 2    30     42.44      0.05123
## 3    40     55.17      0.04468
## 4    50     66.31      0.05029
## 
## 
## Elapsed time : 0.2889349
# Extract the conditional volatility series
conditional_volatility <- sigma(garch_model)

# Plot the conditional volatility
plot(conditional_volatility, type = "l", col = "blue", 
     xlab = "Date", ylab = "Conditional Volatility", 
     main = "Conditional Volatility (GARCH(1,1))")

# GJRGARCH ###########

# Specify the GJR GARCH model without external regressors
spec <- ugarchspec(mean.model = list(armaOrder = c(1, 0)),
                   variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)))

returns_df <- data.frame(Date = data$Date[-1], Returns = log_returns)
# Estimate the model
gjr_garch_model <- ugarchfit(data = returns_df$Returns, spec = spec)

# Print the model summary
show(gjr_garch_model)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001723    0.001281  1.34572 0.178393
## ar1     0.043737    0.044672  0.97908 0.327541
## omega   0.000113    0.000039  2.90256 0.003701
## alpha1  0.285825    0.075806  3.77046 0.000163
## beta1   0.693616    0.059561 11.64554 0.000000
## gamma1 -0.094777    0.085879 -1.10361 0.269763
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001723    0.001252  1.37605 0.168806
## ar1     0.043737    0.044405  0.98496 0.324643
## omega   0.000113    0.000063  1.80367 0.071283
## alpha1  0.285825    0.077750  3.67618 0.000237
## beta1   0.693616    0.081029  8.56012 0.000000
## gamma1 -0.094777    0.095150 -0.99608 0.319211
## 
## LogLikelihood : 1220.442 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -3.9430
## Bayes        -3.8999
## Shibata      -3.9432
## Hannan-Quinn -3.9262
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01619  0.8988
## Lag[2*(p+q)+(p+q)-1][2]   0.15131  0.9996
## Lag[4*(p+q)+(p+q)-1][5]   1.78035  0.7734
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.06384  0.8005
## Lag[2*(p+q)+(p+q)-1][5]   0.91333  0.8790
## Lag[4*(p+q)+(p+q)-1][9]   2.62971  0.8182
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1474 0.500 2.000  0.7010
## ARCH Lag[5]    1.0839 1.440 1.667  0.7082
## ARCH Lag[7]    1.5308 2.315 1.543  0.8152
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.5967
## Individual Statistics:              
## mu     0.22769
## ar1    0.02798
## omega  0.31376
## alpha1 0.31278
## beta1  0.47534
## gamma1 0.44370
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.41711 0.6767    
## Negative Sign Bias 0.24668 0.8052    
## Positive Sign Bias 0.06494 0.9482    
## Joint Effect       0.31421 0.9573    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     29.13      0.06397
## 2    30     37.28      0.13922
## 3    40     47.64      0.16156
## 4    50     63.38      0.08122
## 
## 
## Elapsed time : 0.8259881
# Extract the conditional volatility series
conditional_volatility <- sigma(gjr_garch_model)

# Plot the conditional volatility
plot(conditional_volatility, type = "l", col = "purple", 
     xlab = "Date", ylab = "Conditional Volatility", 
     main = "Conditional Volatility (GJR GARCH)")

# Q9 ########
if (!requireNamespace("knitr", quietly = TRUE)) {
  install.packages("knitr")
}

library(knitr)

# Calculate information criteria for garch_model
info_criteria_garch <- unlist(infocriteria(garch_model))

# Calculate information criteria for gjr_garch_model
info_criteria_gjr_garch <- unlist(infocriteria(gjr_garch_model))

# Create a data frame for the table
info_criteria_table <- data.frame(
  Model = c("GARCH", "GJR-GARCH"),
  Akaike = c(info_criteria_garch[1], info_criteria_gjr_garch[1]),
  Bayes = c(info_criteria_garch[2], info_criteria_gjr_garch[2]),
  Shibata = c(info_criteria_garch[3], info_criteria_gjr_garch[3]),
  HQ = c(info_criteria_garch[4], info_criteria_gjr_garch[4])
)

# Function to determine model type based on criterion name
get_model_type <- function(criterion_values) {
  model_types <- c("GARCH", "GJR")
  selected_model <- model_types[which.min(criterion_values)]
  return(selected_model)
}

# Add a row for the model type
model_types <- apply(info_criteria_table[, -1], 2, get_model_type)
info_criteria_table <- rbind(
  info_criteria_table,
  c("Best Model", model_types)
)

# Print the table
print(kable(info_criteria_table, format = "markdown", digits = 5))
## 
## 
## |Model      |Akaike            |Bayes             |Shibata           |HQ                |
## |:----------|:-----------------|:-----------------|:-----------------|:-----------------|
## |GARCH      |-3.9442934538725  |-3.90839047527239 |-3.94442381253969 |-3.9303336527782  |
## |GJR-GARCH  |-3.94299377098951 |-3.89991019666938 |-3.94318108753535 |-3.92624200967634 |
## |Best Model |GARCH             |GARCH             |GARCH             |GARCH             |