Dataset: Arch Capital Group Ltd. (ACGL) Stock Prices on Kaggle

library(PerformanceAnalytics)
library(quantmod)
library(ggplot2)
library(forecast)
library(tseries)
library(FinTS)
library(rugarch)
library(dplyr)
library(gridExtra)
library(stats)
library(fGarch)
library(xts)
library(evir)

data<-read.csv("acgl.us.csv", header = TRUE, sep = ',')
data <- na.omit(data)
data$Date <- as.Date(data$Date)

closing_prices <- data$Close
closing_xts <- xts(closing_prices, order.by = data$Date)

1 Exploratory Data Analysis

par(mfrow = c(2,2))

hist(closing_prices, main = "Distribution of Closing Prices", col = "steelblue", xlab = "Price")
plot(closing_prices, type='l', main = "Line Plot of the Closing Prices", col = "blue", ylab = "Closing Prices")

acf(closing_prices,  main = "ACF Plot of the Closing Prices")
pacf(closing_prices,  main = "PACF Plot of the Closing Prices")

plot(decompose(ts(closing_prices, frequency = 30)))

summary(closing_prices)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.18   22.64   33.15   41.83   57.55  102.38

2 Assessing Stationarity

Closing Prices have increasing mean and standard deviation. Using log returns to achieve stationarity.

returns <- na.omit(diff(log(closing_prices)))
returns_xts <- xts(returns, order.by = data$Date[-1])
adf_result <- adf.test(returns)

par(mfrow = c(2,2))

hist(returns, main = "Distribution of Log Returns", col = "dodgerblue3", xlab = "Returns")
plot(returns,  main = "Log Returns Over Time", col = "darkred", ylab = "Returns")

acf(returns,  main = "ACF Plot of the Returns")
pacf(returns,  main = "PACF Plot of the Returns", xlab = "Value")

adf_result
## 
##  Augmented Dickey-Fuller Test
## 
## data:  returns
## Dickey-Fuller = -16.56, Lag order = 14, p-value = 0.01
## alternative hypothesis: stationary

3 Fitting a Seasonal ARIMA Model

Identifying an appropriate SARIMA model and validating its fit by checking the residuals.

best_fit <- auto.arima(returns, seasonal = TRUE)

par(mfrow = c(1,2))
plot(residuals(best_fit), type = "l", main = "Line Plot of SARIMA Residuals", ylab = "Residuals", xlab = "Time", col = "blue")
abline(h = 0, col = "red", lty = 2)

acf(residuals(best_fit), main = "ACF of ARIMA Residuals")

summary(best_fit)
## Series: returns 
## ARIMA(4,0,3) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3     ar4      ma1     ma2      ma3   mean
##       0.1231  -0.2504  0.8296  0.0623  -0.2360  0.2259  -0.8502  6e-04
## s.e.  0.0689   0.0668  0.0455  0.0229   0.0663  0.0736   0.0467  1e-04
## 
## sigma^2 = 0.0001751:  log likelihood = 9303.4
## AIC=-18588.8   AICc=-18588.74   BIC=-18534.16
## 
## Training set error measures:
##                         ME       RMSE         MAE MPE MAPE      MASE
## Training set -3.372768e-06 0.01321677 0.008973802 NaN  Inf 0.6694041
##                     ACF1
## Training set 0.001227962

4 Modeling Volatility Clustering

arch_test <- ArchTest(residuals(best_fit), lags = 12)
arch_test
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  residuals(best_fit)
## Chi-squared = 705.02, df = 12, p-value < 2.2e-16

Using the GARCH(1,1) Model with Student-t Distribution to model the ARCH effects.

Validating the model by checking squared ACF for remaining ARCH effects and plotting standardized residuals.

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

fit <- ugarchfit(data = returns, spec = spec)

par(mfrow = c(1,2))

resid_std <- residuals(fit, standardize = TRUE)
plot(resid_std, type = "l", main = "Standardized Residuals from GARCH(1,1)",
     ylab = "Residuals", xlab = "Time", col = "darkblue")
acf(resid_std^2, main = "ACF of Squared Standardized Residuals")

5 Extreme Value Analysis

Forecasting volatility one day ahead and computing Value-at-Risk (VaR) and Expected Shortfall (ES).

forecast <- ugarchforecast(fit, n.ahead = 1)
sigma <- sigma(forecast)

probs <- c(0.95, 0.99, 0.999)
df <- coef(fit)["shape"]
mu <- 0

VaR <- -qt(probs, df = df) * sigma + mu
ES  <- -sigma * (dt(qt(probs, df), df) / (1 - probs)) *
  ((df + qt(probs, df)^2) / (df - 1)) + mu

VaR <- VaR * 100
ES <- ES * 100

risk_table <- data.frame(
  Confidence = probs,
  VaR = round(VaR, 5),
  ES = round(ES, 5)
)

ggplot(risk_table, aes(x = as.factor(Confidence))) +
  geom_bar(aes(y = VaR, fill = "VaR"), stat = "identity", position = "dodge") +
  geom_point(aes(y = ES, color = "ES"), size = 3) +
  scale_fill_manual("", values = c("VaR" = "dodgerblue")) +
  scale_color_manual("", values = c("ES" = "darkred")) +
  labs(title = "Value-at-Risk (VaR) and Expected Shortfall (ES)",
    x = "Confidence Level", y = "Potential Loss of Portfolio Value (%)") +
  theme_minimal(base_size = 13)

risk_table
##   Confidence      VaR       ES
## 1      0.950 -2.00136 -2.80279
## 2      0.990 -3.25154 -4.19019
## 3      0.999 -5.43054 -6.72264