Pre-Processing

#Set Working direc to call file
setwd("C:\\Users\\javen\\OneDrive\\Desktop\\School\\RM 708")

#read data
library(readxl)
data <- read_excel("C:\\Users\\javen\\OneDrive\\Desktop\\School\\RM 708\\GOOG2.xlsx", col_names=TRUE)

#Make sure data loaded in appropriately 
head(data, 10)
## # A tibble: 10 x 7
##    Date                 Open  High   Low Close Adj_Close      Volume
##    <dttm>              <dbl> <dbl> <dbl> <dbl>     <dbl>       <dbl>
##  1 2004-10-01 00:00:00  3.26  4.98  3.21  4.75      4.75 10321258145
##  2 2004-11-01 00:00:00  4.82  5.02  4.02  4.53      4.53 11145383635
##  3 2004-12-01 00:00:00  4.53  4.98  4.20  4.80      4.80  5835483147
##  4 2005-01-01 00:00:00  4.92  5.11  4.39  4.87      4.87  8237356337
##  5 2005-02-01 00:00:00  4.84  5.40  4.51  4.68      4.68 13073210618
##  6 2005-03-01 00:00:00  4.71  4.73  4.30  4.50      4.50  6903740285
##  7 2005-04-01 00:00:00  4.53  5.60  4.48  5.48      5.48  9101631755
##  8 2005-05-01 00:00:00  5.53  6.93  5.48  6.91      6.91 10188799510
##  9 2005-06-01 00:00:00  7.05  7.70  6.66  7.33      7.33 15618390955
## 10 2005-07-01 00:00:00  7.35  7.92  7.15  7.17      7.17  8228691980
# Calculate returns
library(dplyr)

# Convert date to Date type
data$Date <- as.Date(data$Date, format = "%m/%d/%Y")

# Arrange data by date and calculate returns
data <- data %>%
  arrange(Date) %>%
  mutate(Return = 100 * (Close - lag(Close)) / lag(Close))

# Arrange data by date and calculate returns
data <- data %>%
  arrange(Date) %>%
  mutate(Adj_Return = 100 * (Adj_Close - lag(Adj_Close)) / lag(Adj_Close))

# Append Log(Adj Close)
data <- data %>%
  mutate(log_Adj_Close = log(Adj_Close))

# Remove NA values
data <- na.omit(data)

#Convert to time-series
library(TSA)
GOOG2 <- ts(data$Return, start=c(2004, 10), frequency=12)

#Make sure Returns were calculated appropriately 
head(GOOG2, 10)
##  [1] -4.5426027  5.9402172  1.4679225 -3.9004252 -3.9789415 21.8769198
##  [7] 26.0318075  6.0879267 -2.1723549 -0.6116251
\[\text{}\]

\[\text{I. Consider the following two MA models for a sequence } y_t\] \[\text{(i) } y_t = u_t + \theta_1u_{t-1}\] \[\text{(ii) }y_t = u_t + \theta_1u_{t-1} + \theta_2u_{t-2}\]

\[\text{Where } u_t \text{ is a sequence of iid random variables with } E(u_t)=0 \text{ and } Var(u_t) = \sigma^2, \text{at any time t}.\]

\[\text{}\]

\(\text{a) Calculate the mean and variance of } y_t \text{ from model (i).}\) \[\text{}\] \[\text{Estimate } y_t = u_t + \theta_1u_{t-1}\]

MA1 <- arima(GOOG2, order = c(0,0,1))
MA1fitted_values <- fitted(MA1)
MA1mean_estimated <- round(mean(MA1fitted_values),5)
MA1variance_estimated <- round(var(MA1fitted_values),5)
## Estimated Mean: 1.77671
## Estimated Variance: 2e-05

\(\text{ b) Derive the autocorrelation function for the process in model (i).}\)

\[\text{}\]

\[y_t = u_t + \theta_1 u_{t-1}\]

\[\text{Let the autocorrelation function at lag } k \text { be } \rho(k)\].

\[\rho(k) = \frac{Cov(y_t, y_{t-k})}{Var(y_t)}\]

\[\text{When } k = 0\]

\[\rho(0) = \frac{Cov(y_t, y_t)}{Var(y_t)} = 1\]

\[\text{When } k = 1\]

\[\rho(1) = \frac{Cov(y_t, y_{t-1})}{Var(y_t)} = \frac{E[(u_t + \theta_1u_{t-1})(u_{t-1} + \theta_1u_{t-2})]}{(1 + \theta_1^2)\sigma^2} = \frac{\theta_1\sigma^2}{(1 + \theta_1^2)\sigma^2} = \frac{\theta_1}{1 + \theta_1^2}\]

\[\text{When } k > 1\]

\[\rho(k) = 0 \text{, as there is no correlation between }u_t \text{ and } u_{t-h} \text{ for } h > 1 \text{ in a MA(1) process.}\]

\[\text{So, the ACF of a MA(1) process is 1 at lag 0, and } \frac{\theta_1}{1 + \theta_1^2} \text{ at lag 1, and 0 otherwise.}\] >\[\text{Estimate on } y_t = u_t + \theta_1u_{t-1} \text{ residuals}\]

acf(MA1$residuals, lag.max = 48)

\(\text{c) Calculate the mean and variance of} y_t \text{from model (ii).}\) \[\text{}\] \[\text{Estimate } y_t = u_t + \theta_1u_{t-1} + \theta_2u_{t-2}\]

MA2 = arima(GOOG2, order=c(0,0,2))
MA2fitted_values <- fitted(MA2)
MA2mean_estimated <- round(mean(MA2fitted_values),5)
MA2variance_estimated <- round(var(MA2fitted_values),5)
## Estimated Mean: 1.77767
## Estimated Variance: 0.26942

\(\text{d) Derive the autocorrelation function for the process in model (ii).}\)

\[\text{}\]

\[y_t = u_t + \theta_1u_{t−1} + \theta_2u_{t−2}\].

\[\text{Similarly, the ACF for lag } k \text{here is given by:}\]

\[\rho(k) = \frac{Cov(y_t, y_{t-k})}{Var(y_t)}\]

\[ \text{When }k = 0\]

\[\rho(0) = \frac{Cov(y_t, y_t)}{Var(y_t)} = 1\]

\[ \text{When } k = 1\]

\[\rho(1) = \frac{Cov(y_t, y_{t-1})}{Var(y_t)} = \frac{E[(u_t + \theta_1u_{t-1} + \theta_2u_{t-2})(u_{t-1} + \theta_1u_{t-2} + \theta_2u_{t-3})]}{(1 + \theta_1^2 + \theta_2^2)\sigma^2} = \frac{(\theta_1 + \theta_2)\sigma^2}{(1 + \theta_1^2 + \theta_2^2)\sigma^2} = \frac{(\theta_1 + \theta_2)}{1 + \theta_1^2 + \theta_2^2}\] \[\text{When } k = 2\]

\[\rho(2) = \frac{Cov(y_t, y_{t-2})}{Var(y_t)}= \frac{E[(u_t + \theta_1u_{t-1} + \theta_2u_{t-2})(u_{t-2} + \theta_1u_{t-3} + \theta_2u_{t-4})]}{(1 + \theta_1^2 + \theta_2^2)\sigma^2} = \frac{\theta_2\sigma^2}{(1 + \theta_1^2 + \theta_2^2)\sigma^2} = \frac{\theta_2}{1 + \theta_1^2 + \theta_2^2}\]

\[\text{When } k > 2\]

\[\rho(k) = 0 \text{, because there's no correlation between} u_t \text{ and } u_{t-h} \text{ for } h > 2 \text{ in a MA(2) process.}\]

\[\text{So, the ACF of a MA(2) process is 1 at lag 0, } \frac{(\theta_1 + \theta_2)}{1 + \theta_1^2 + \theta_2^2} \text{at lag 1, } \frac{\theta_2}{1 + \theta_1^2 + \theta_2^2} \text{at lag 2, and 0 otherwise.}\] \[\text{}\]

\[\text{Estimate on } y_t = u_t + \theta_1u_{t-1} + \theta_2u_{t-2} \text{ residuals}\]

acf(MA2$residuals, lag.max=48)

\(\text{(e) Using the model (ii), calculate all the autocorrelation coefficients if } \theta_1=-0.4 , \text{ and } \theta_2=0.3.\) \[\rho(0) = 1\] \[\rho(1) = \frac{-.4+.3}{1+.4^2+.03^2} \frac{-.1}{1+.16+.09} = \frac{-.1}{1.25} = -.08\] \[\rho(2) = \frac{.3}{1.25} = .24\] \[\rho(k>2) = 0\]

\[\text{II. In this exercise you will use time series analysis techniques to model and forecast monthly GOOGLE price data,}\]

\[\text{using the ‘Adjusted Closed’ price column of the attached GOOG file. }\]

\[\text{}\]

\(\text{a) Load the GOOG data file attached to this question file and examine its structure,}\) \(\text{on the basis of the information you will generate after having done the following:}\) \[\text{}\] \[\text{(i) Plot the log(price)}\]

#I appended Log(Adj_Close) when loading the data - see top of file
GOOG2_LOG_Adj_Close <- ts(data$log_Adj_Close, start=c(2004, 10), frequency=12)
plot.ts(GOOG2_LOG_Adj_Close, 
        xlab = "Date", 
        ylab = "Log of Adj Close", 
        main = "Log of Google's Adjusted Close Price Over Time", 
        col = "blue")

\[\text{(ii) Plot the ACF as well as the PACF of log(price)}\]

acf(GOOG2_LOG_Adj_Close, lag.max=150)

\[\text{The observed ACF seems to dampen around lag 4}\] \[\text{}\]

pacf(GOOG2_LOG_Adj_Close, lag.max=150)

\[\text{The observed PACF seems to dampen around lag 2}\] \[\text{}\]

\[\text{(iii) Comment and specify what order of ARIMA model you would recommend.}\] \[\text{}\] \[\text{Based on the observed behavior of the ACF and PACF plots,}\] \[\text{the Log(Adj Closing Price) residuals follow AR order = 4 and MA order = 2.}\] \[\text{Thus, I reccomend ARIMA parameters of (4,0,2) }\]

\[\text{b) Fit and examine the ARIMA model recommended in question a). }\]

ARIMA_4_0_2 <- arima(GOOG2_LOG_Adj_Close, order = c(4,0,2))
ARIMA_4_0_2_fitted_values <- fitted(ARIMA_4_0_2)
summary(ARIMA_4_0_2_fitted_values)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.506   2.551   3.268   3.274   4.006   5.009
## Estimated Variance: 0.7926

\[\text{c) Fit and examine the ARIMA model recommended in question a). }\]

# Extract the residuals
residuals <- residuals(ARIMA_4_0_2)
# Perform a Ljung-Box test to check for autocorrelation in the residuals
Box.test(residuals, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 0.26143, df = 1, p-value = 0.6091

\[\text{}\] \[\text{The high p-value observed here is strong evidence}\] \[\text{that the residuals can be considered indpendent.}\]

\[\text{}\]

# Plot the residuals
plot(residuals, main="Residuals")

\[\text{The randomness with relatively constant variance observed here }\]

\[\text{is furthermore soldifying evidencethat the residuals can be considered indpendent.}\] \[\text{}\]

# Plot the ACF of the residuals
acf(residuals, main="ACF of Residuals")

# Plot the PACF of the residuals
pacf(residuals, main="PACF of Residuals")

\[\text{In both the ACF and PACF, we see all bars within the confidence band.}\]

\[\text{This strengthens the hypothesis that the residuals can be considered indpendent.}\] \[\text{And this furthermore suggests that the correct AR and MA orders were selected.}\] \[\text{}\]

# QQ plot of the residuals to check normality
qqnorm(residuals, main="QQ Plot of Residuals")
qqline(residuals)

\[\text{Here we see the residuals follow a 45-degree reference line quite closely.}\] \[\text{This purports well-behaved, normal residuals.}\]

\[\text{d) Make a 24 months ahead forecast of log(price) series, and if possible, comment on your projections. }\]

library(forecast)
# Make a 24-month ahead forecast
forecast <- predict(ARIMA_4_0_2, n.ahead = 24)

# The forecast object contains predicted values and standard errors
predicted_values <- forecast$pred
se <- forecast$se

# Calculate the 95% prediction intervals
upper_bound <- predicted_values + 2*se
lower_bound <- predicted_values - 2*se

# Increment the period by 1 to get the start time of the forecasted values
# Get the end time of the original series
end_time <- end(GOOG2_LOG_Adj_Close)
start_time <- end_time
start_time[2] <- start_time[2] + 1
if (start_time[2] > 12) { # If the period exceeds 12, increment the year by 1 and set the period to 1
  start_time[1] <- start_time[1] + 1
  start_time[2] <- 1
}
# Create a time series for the forecasted values with the correct time points
predicted_values_ts <- ts(forecast$pred, start = start_time, frequency = 12)

# Create a time series for the upper and lower bounds of the prediction intervals
upper_bound_ts <- ts(upper_bound, start = start_time, frequency = 12)
lower_bound_ts <- ts(lower_bound, start = start_time, frequency = 12)

# Plot the original series, forecast and prediction intervals
plot(GOOG2_LOG_Adj_Close, xlim = c(start(GOOG2_LOG_Adj_Close)[1], end(predicted_values_ts)[1]), main = "Original Series, Forecast and Prediction Intervals")
lines(predicted_values_ts, col = "blue")
lines(upper_bound_ts, col = "red", lty = 2)
lines(lower_bound_ts, col = "red", lty = 2)

# Get the time point of the last point of the original series
end_time <- as.numeric(end(GOOG2_LOG_Adj_Close)[1]) + (end(GOOG2_LOG_Adj_Close)[2] - 1)/12

# Get the time point of the first point of the forecasted values
start_time <- as.numeric(start(predicted_values_ts)[1]) + (start(predicted_values_ts)[2] - 1)/12

# Add a line segment that connects the last point of the original series to the first point of the forecasted values
lines(c(end_time, start_time), c(tail(GOOG2_LOG_Adj_Close, n = 1), forecast$pred[1]), col = "blue")

\[\text{It is observed that the 95% confidence interval (dashed-arch in red)}\]

\[\text{is within reasonable distance of the relative min and max for the past 5 years.}\] \[\text{The prediction line in blue also seens to follow the general shape of the curve}\]

\[\text{and does not appear to be biased towards any noise.}\]

\[\text{The predictions therfore seem reasonable.}\]

\[\text{}\]

\[\text{III. In this exercise you will analyze the volatility of monthly GOOGLE returns, or log(returns).}\]

\[\text{a) Load the GOOG data file attached to this question file and examine its structure,}\] \[\text{on the basis of the information you will generate after having done the following:}\]

\[\text{(i) Use the GOOGLE adjusted closing prices to calculate google monthly returns.}\] \[\text{This was done at the top of the file in the pre-processing stage.}\] \[\text{Below is a sample.}\]

GOOG2_Adj_Return = ts(data$Adj_Return, start=c(2004, 10), frequency=12)
head(GOOG2_Adj_Return)
##            Jan       Feb       Mar Apr May Jun Jul Aug Sep       Oct       Nov
## 2004                                                       -4.542603  5.940217
## 2005 -3.900425 -3.978942 21.876920                                            
##            Dec
## 2004  1.467923
## 2005

\[\text{(ii) Plot the series of returns over the time.}\]

plot.ts(GOOG2_Adj_Return, 
        xlab = "Date", 
        ylab = "Adj Return", 
        main = "Google's Adjusted Return Over Time", 
        col = "blue")

\[\text{(iii) Plot the ACF and the PACF of the returns, as well as of the absolute value of the returns.}\]

#Full ACF view on Adj_Return
acf(GOOG2_Adj_Return, lag.max = 300)

#Relevant ACF Lags on Adj_Returns
acf(GOOG2_Adj_Return, lag.max = 60)

#Full PACF view on Adj_Return
pacf(GOOG2_Adj_Return, lag.max = 300)

#Relevant PACF Lags on Adj_Returns
pacf(GOOG2_Adj_Return, lag.max = 60)

#Full ACF view on ABS(Adj_Return)
acf(abs(GOOG2_Adj_Return), lag.max = 300)

#Relevant ACF Lags on Adj_Returns
acf(abs(GOOG2_Adj_Return), lag.max = 60)

#Full PACF view on ABS(Adj_Return)
pacf(abs(GOOG2_Adj_Return), lag.max = 300)

#Relevant PACF Lags on Adj_Returns
pacf(abs(GOOG2_Adj_Return), lag.max = 60)

\[\text{(iv) Plot the ACF and the PACF of the squared returns.}\]

#Full ACF view on Squared Adj_Returns
acf(GOOG2_Adj_Return^2, lag.max = 300)

#Relevant ACF Lags on Squared Adj_Returns
acf(GOOG2_Adj_Return^2, lag.max = 60)

#Full PACF View on Squared Adj_Returns
pacf(GOOG2_Adj_Return^2, lag.max = 300)

#Relevant PACF Lags on Squared Adj_Returns
pacf(GOOG2_Adj_Return^2, lag.max = 60)

\[\text{(v) Compute the mean of the returns,}\] \[\text{and test whether the mean return is equal to zero or greater,}\] \[\text{with a 5% significance level.}\]

# Compute the mean of the returns
mean_return = mean(GOOG2_Adj_Return)

cat("Mean:",mean_return)
## Mean: 1.776732
# Perform a one-sample t-test to test whether the mean return is equal to zero
t_test = t.test(GOOG2_Adj_Return, mu = 0)

# Check the p-value
if(t_test$p.value < 0.05) {
  print("There is sufficient evidence to reject the null hypothesis that the mean return is equal to zero.")
} else {
  print("There is not sufficient evidence to reject the null hypothesis that the mean return is equal to zero.")
}
## [1] "There is sufficient evidence to reject the null hypothesis that the mean return is equal to zero."

\[\text{(b) Test for the ARCH effect in the mean-adjusted squared returns}\]

# Perform a one-sample one-sided t-test to test whether the mean return is greater than zero
t_test_greater = t.test(GOOG2_Adj_Return, mu = 0, alternative = "greater")

# Check the p-value
if(t_test_greater$p.value < 0.05) {
  print("There is sufficient evidence to reject the null hypothesis that the mean return is equal to or less than zero.")
} else {
  print("There is not sufficient evidence to reject the null hypothesis that the mean return is equal to or less than zero.")
}
## [1] "There is sufficient evidence to reject the null hypothesis that the mean return is equal to or less than zero."
# Calculate the residuals (the difference between actual values and the mean)
residuals = GOOG2_Adj_Return - mean(GOOG2_Adj_Return)

# Square the residuals
squared_residuals = residuals^2

# Conduct the Ljung-Box test on the squared residuals
Box.test(squared_residuals, lag = 12, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  squared_residuals
## X-squared = 27.398, df = 12, p-value = 0.00677

\[\text{Here a p-value of less than 1% is observed.}\] \[\text{It is therefore highly likely that}\] \[\text{ARCH effects are present in the mean-adjusted squared returns.}\]

\[\text{(c) Fit and examine the GARCH model that you believe is appropriate}\] \[\text{for return or log(return) series.}\]

library(rugarch)

# Specify the GARCH model
spec = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,2)))

# Fit the GARCH model
fit = ugarchfit(spec, data = GOOG2_Adj_Return)

# Print the model
print(fit)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,2)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error     t value Pr(>|t|)
## mu       1.80475    0.043892    41.11744 0.000000
## ar1      0.93242    0.019706    47.31759 0.000000
## ma1     -0.99605    0.000764 -1303.46961 0.000000
## omega    2.38880    2.389084     0.99988 0.317368
## alpha1   0.10144    0.052272     1.94054 0.052314
## beta1    0.13937    0.149022     0.93524 0.349663
## beta2    0.72121    0.145239     4.96569 0.000001
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu       1.80475    0.042077    42.8919 0.000000
## ar1      0.93242    0.018797    49.6052 0.000000
## ma1     -0.99605    0.000910 -1094.9574 0.000000
## omega    2.38880    2.064394     1.1571 0.247214
## alpha1   0.10144    0.039980     2.5371 0.011176
## beta1    0.13937    0.091668     1.5204 0.128409
## beta2    0.72121    0.079565     9.0644 0.000000
## 
## LogLikelihood : -775.6492 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       7.0193
## Bayes        7.1262
## Shibata      7.0174
## Hannan-Quinn 7.0625
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0003437  0.9852
## Lag[2*(p+q)+(p+q)-1][5] 0.3520623  1.0000
## Lag[4*(p+q)+(p+q)-1][9] 1.2201100  0.9992
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                     0.05915  0.8078
## Lag[2*(p+q)+(p+q)-1][8]    3.97670  0.5103
## Lag[4*(p+q)+(p+q)-1][14]   6.95158  0.5121
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]     1.844 0.500 2.000  0.1745
## ARCH Lag[6]     4.782 1.461 1.711  0.1267
## ARCH Lag[8]     5.205 2.368 1.583  0.2278
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.7337
## Individual Statistics:              
## mu     0.07256
## ar1    0.33026
## ma1    0.13023
## omega  0.10328
## alpha1 0.14018
## beta1  0.12679
## beta2  0.12293
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.5325 0.5949    
## Negative Sign Bias  1.0361 0.3013    
## Positive Sign Bias  0.2047 0.8380    
## Joint Effect        1.1201 0.7722    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     20.32       0.3757
## 2    30     28.17       0.5091
## 3    40     35.83       0.6151
## 4    50     47.18       0.5472
## 
## 
## Elapsed time : 0.115566

\[\text{The } \mu, \text{ 'ar1', and 'ma1' parameters are statistically significant,}\] \[\text{as their p-values are 0.000000 which is less than any common level of significance (0.05, 0.01, 0.001).}\]

\[\text{The } \alpha_1 \text{ and } \beta_2 \text{ parameters are also statistically significant,}\] \[\text{but } \omega \text{ and } \beta_1 \text{ parameters are not significant according to the regular standard errors.}\]

\[\text{When considering the robust standard errors,}\] \[\omega \text{ and } \beta_1 \text{ remain insignificant, but } \alpha_1 \text{ becomes statistically significant - however.}\]

\[\text{The log-likelihood of the model is -775.6492.}\] \[\text{Lower absolute values are generally better - but this is not bad,}\] \[\text{and this metric is most useful when comparing multiple models anyhow.}\]

\[\text{The Information criteria (Akaike, Bayes, Shibata, Hannan-Quinn)}\] \[\text{are also used for model comparison - the lower these values are,}\] \[\text{the better the model is considered in terms of its trade-off between goodness of fit and model complexity.}\] \[\text{With approximate 7s across the board, I don't see many models scoring much lower than this.}\]

\[\text{The Weighted Ljung-Box Test and Weighted ARCH LM Tests have high p-values,}\] \[\text{suggesting that there is no significant serial correlation in the residuals or their squares,}\] \[\text{which means the model is a good fit.}\]

\[\text{The Nyblom stability test suggests that the parameters of the model are stable over time.}\]

\[\text{The Sign Bias Test suggests there's no significant sign bias in residuals.}\]

\[\text{The Adjusted Pearson Goodness-of-Fit Test has high p-values suggesting that the model fits the data well.}\]

\[\text{(d) Examine the residuals as well as the volatility}\] \[\text{of the GARCH model estimated in c).}\]

# Residuals Plot
plot(residuals(fit), main="Residuals of the GARCH model", ylab="Residuals")

# Residuals should look like white noise, that is, they should be randomly dispersed around the zero line. 
# Observed patterns in the residuals could suggest that there is more information in the data that the model is not capturing.

# Standardized Residuals Over Time
plot(residuals(fit, standardize=T), main="Standardized Residuals of the GARCH model", ylab="Standardized Residuals")

\[\text{}\]

\[\text{Both the pure and standardized residuals appear as white noise,}\]

\[\text{suggesting that there isn't any explanatory information that the model is not capturing.}\]

\[\text{}\]

# Density Plot of Residuals
plot(density(residuals(fit)), main="Density of Residuals", xlab="Residuals")
lines(density(rnorm(1000, mean(residuals(fit)), sd(residuals(fit)))), lty="dotted", col="red")

# This plot should look approximately like a bell curve if the residuals are normally distributed. 
# If the distribution is significantly skewed or has heavy tails, it could suggest that the model's assumptions are not being met.

# Density Plot of Standardized Residuals
plot(density(residuals(fit, standardize=T)), main="Density of Standardized Residuals", xlab="Standardized Residuals")
lines(density(rnorm(1000, mean(residuals(fit, standardize=T)), sd(residuals(fit, standardize=T)))), lty="dotted", col="red")

\[\text{Both residual density distributions (black line) follows}\] \[\text{ a normal distribution with the same mean and standard deviation as the residuals(red line) quite closely.}\] \[\text{Furthermore, the distribution does not appear to be significantly skewed,}\] \[\text{and there are no especially heavy tails observed,}\]

\[\text{suggesting that GARCH(1,2) assumptions are met.}\] \[\text{}\]

# Q-Q Plot of Residuals
qqnorm(residuals(fit), main="Q-Q Plot of Residuals")
qqline(residuals(fit))

# In a Q-Q plot, the quantiles of the residuals are plotted against the quantiles of a normal distribution. 
# If the residuals are normally distributed, the points should fall approximately along the reference line.

# Q-Q Plot of Standardized Residuals
qqnorm(residuals(fit, standardize=T), main="Q-Q Plot of Standardized Residuals")
qqline(residuals(fit, standardize=T))

\[\text{The vast majority of residuals observed in both Q-Q plots follow a 45-degree refrence line quite closely.}\] \[\text{This furthermore suggests that the model residuals are normally distributed.}\]

\[\text{}\]

# Volatility Plot
plot(sigma(fit), main="Volatility of the GARCH model", ylab="Volatility")

# This plot shows the estimated volatility over time. In many financial time series, volatility clusters in certain periods (volatility clustering), which is something GARCH models can capture.

\[\text{Here a fairly typical volitality graph is observed given the data at hand is stock return data.}\]