QUESTION ONE.

#Loading the required packages
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## 
## The following object is masked from 'package:graphics':
## 
##     legend
#Importing the data set
stock.price <-read.csv("C:/Users/Baha/Downloads/HistoricalPrices.csv")
#Selecting the Close stock price column.
stock.price <-stock.price[,-c(2,3,4)]
#Calculate log returns
log_returns <- diff(log(stock.price$Close))
#Create a time series plot for the stock
plot.ts(log_returns, main = "Time series plot of Stock S&P500")

Conclusion. The time series plot of log returns of the stock is stationary. This implies that the mean of log returns is not different from zero.

#summary statistics of the stock log returns
stacs <- data.frame(
  "Close_Stock_price" = c(mean(log_returns),sd(log_returns),skewness(log_returns),kurtosis(log_returns, type = "excess"),
         min(log_returns),max(log_returns)))
row.names(stacs) <-c("mean","standard deviation","skewness","kurtosis","minimum","maximum")
print(stacs)
##                    Close_Stock_price
## mean                    0.0003108078
## standard deviation      0.0120831263
## skewness               -0.5095683866
## kurtosis               12.7629190291
## minimum                -0.1276521412
## maximum                 0.1095719593

Conclusion. The mean of log returns is zero, implying that the log returns are stationary. Further, from the statistics, the skewness is 0.51(2.dp), and for Normal Distribution Skewness is always zero. This is evident that the log returns are not normally distributed. Also, there is presence of heavy tails, this is concluded from the calculated excess kurtosis. To deal with such data, one may use Frechet Distribution or Student-t distribution.

#Testing for the presence of the mean effect in the log returns
#perform a two-sided t-test for log returns of the stock
t.test <-t.test(log_returns, mu = 0)
#extract the p-values for each t.test
p.value <-c(t.test$p.value)
print(p.value)
## [1] 0.06438895

Conclusion. Since the p-value is greater than 0.05, we fail to reject the null hypothesis at the 5% significance level. This means that there is not enough evidence to conclude that the mean log return of the close stock is statistically different from zero. However, it is important to note that the absence of evidence for a difference does not necessarily imply evidence of absence of a difference. Therefore, further research may be needed to draw a more definitive conclusion.

library(tseries)
test_result <-jarque.bera.test(stock.price$Close)
test_statistic <- test_result$statistic
p.value <-test_result$p.value
#print the results
cat("Jarque-Bera test statistic:", test_statistic, "\n")
## Jarque-Bera test statistic: 856.8553
cat("p.value:", p.value, "\n")
## p.value: 0
#checking the significance at 5% level
if(p.value<0.05){
  cat("Reject the null hypothesis: the log returns are not normally distributed\n")
} else{
  cat("Do not reject the null hypothesis: the log returns are normally distributed\n")
}
## Reject the null hypothesis: the log returns are not normally distributed

Conclusion. This further proves that the log returns are not normally distributed using the Jarque Bera test. The p-value is less than 5% level of significance and there is sufficient evidence to reject the null hypothesis that says the log returns are normally distributed. Thus, the log returns are not normally distributed as proven earlier using the summary statistics, the skewness was 0.51, not different from zero.

#Testing for the presence of ARCH effects.
#We use the ARCH-LM test
library(forecast)
arima <-Arima(log_returns, order = c(2,0,2))
residuals_stock <-residuals(arima)
library(FinTS)
## 
## Attaching package: 'FinTS'
## The following object is masked from 'package:forecast':
## 
##     Acf
arch.effect <-ArchTest(residuals_stock);arch.effect
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  residuals_stock
## Chi-squared = 1577.2, df = 12, p-value < 2.2e-16

Conclusion. From the test, we reject the null hypothesis since the p-value of the log returns is less than 5% level of significance, hence the log returns have ARCH EFFECTS.

#Testing for heavy tails, we perform the Anderson Darling Test
library(nortest)
ad.test(residuals_stock)
## 
##  Anderson-Darling normality test
## 
## data:  residuals_stock
## A = 116.16, p-value < 2.2e-16

Conclusion. Testing at 5% level of significance with the hypotheses; Null hypothesis the residuals follow a normal distribution against the alternative hypothesis the residuals do not follow a normal distribution, we reject the null hypothesis since the p-value is less than 5% level of significance. Thus the residuals are not normally distributed.

#Testing for serial correlation of the log returns of stock
#We perform the Ljung Box Test on the log returns
Box_test <-Box.test(log_returns, lag = 10, type = "Ljung-Box")
print(Box_test)
## 
##  Box-Ljung test
## 
## data:  log_returns
## X-squared = 123.2, df = 10, p-value < 2.2e-16

Conclusion. Testing at 5% level of significance and with the hypotheses such that the null hypothesis is that there is no serial correlation, we reject the null hypothesis since the p-value is less than 5% level of significance and thus the log returns are serially correlated.

#Fitting the ARCH(1) model
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:purrr':
## 
##     reduce
## The following object is masked from 'package:stats':
## 
##     sigma
arch_model <-ugarchspec(variance.model = list(model = "sGARCH", garchOrder =c(1,0)),
                        mean.model = list(armaOrder = c(1,0)), distribution = "norm")
fit_arch <-ugarchfit(arch_model, data = residuals_stock)
print(fit_arch)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,0)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000391    0.000124   3.1512 0.001626
## ar1    -0.117389    0.014351  -8.1798 0.000000
## omega   0.000078    0.000002  34.4455 0.000000
## alpha1  0.534749    0.034429  15.5321 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000391    0.000179   2.1820 0.029110
## ar1    -0.117389    0.051512  -2.2789 0.022675
## omega   0.000078    0.000006  12.6201 0.000000
## alpha1  0.534749    0.107966   4.9529 0.000001
## 
## LogLikelihood : 16061.42 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2094
## Bayes        -6.2043
## Shibata      -6.2094
## Hannan-Quinn -6.2076
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      58.49 2.043e-14
## Lag[2*(p+q)+(p+q)-1][2]     58.86 0.000e+00
## Lag[4*(p+q)+(p+q)-1][5]     64.25 0.000e+00
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      11.13 0.0008489
## Lag[2*(p+q)+(p+q)-1][2]     88.77 0.0000000
## Lag[4*(p+q)+(p+q)-1][5]    245.62 0.0000000
## d.o.f=1
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[2]     155.1 0.500 2.000       0
## ARCH Lag[4]     272.0 1.397 1.611       0
## ARCH Lag[6]     357.0 2.222 1.500       0
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.6911
## Individual Statistics:             
## mu     0.8456
## ar1    0.4266
## omega  2.4261
## alpha1 0.7863
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias           2.4915 0.01275  **
## Negative Sign Bias  0.5824 0.56031    
## Positive Sign Bias  0.3235 0.74634    
## Joint Effect       10.4702 0.01496  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     529.9   2.869e-100
## 2    30     567.6   3.325e-101
## 3    40     599.0   2.125e-101
## 4    50     622.2   2.604e-100
## 
## 
## Elapsed time : 1.440901

Conclusion: The estimated value of alpha 1 indicates the impact of past shocks on volatility of the stock. A higher alpha1 value suggests higher volatility persistence of the stock, while a lower alpha1 value suggests lower volatility persistence of the stock.

#Fitting a GARCH(1,1) model with the normal distribution
#Defining the model specifications
model <-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                   mean.model = list(armaOrder = c(1,0)))
#Fitting the model
model_fit <-ugarchfit(model, data = residuals_stock)
print(model_fit)
## 
## *---------------------------------*
## *          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.000407    0.000113    3.591 0.000329
## ar1     0.056662    0.015122    3.747 0.000179
## omega   0.000002    0.000001    3.273 0.001064
## alpha1  0.125049    0.010011   12.491 0.000000
## beta1   0.854759    0.010682   80.020 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000407    0.000099  4.10461 0.000040
## ar1     0.056662    0.012982  4.36462 0.000013
## omega   0.000002    0.000004  0.57453 0.565608
## alpha1  0.125049    0.026298  4.75514 0.000002
## beta1   0.854759    0.041726 20.48518 0.000000
## 
## LogLikelihood : 16926.55 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.5435
## Bayes        -6.5372
## Shibata      -6.5435
## Hannan-Quinn -6.5413
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.074 0.30001
## Lag[2*(p+q)+(p+q)-1][2]     1.391 0.49944
## Lag[4*(p+q)+(p+q)-1][5]     5.455 0.07361
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.9628  0.3265
## Lag[2*(p+q)+(p+q)-1][5]    4.3275  0.2159
## Lag[4*(p+q)+(p+q)-1][9]    5.4697  0.3638
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3] 0.0001959 0.500 2.000  0.9888
## ARCH Lag[5] 1.8508640 1.440 1.667  0.5051
## ARCH Lag[7] 2.2230946 2.315 1.543  0.6701
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  18.6444
## Individual Statistics:             
## mu     0.1493
## ar1    0.1258
## omega  2.4718
## alpha1 0.1856
## beta1  0.3225
## 
## 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           3.4509 5.633e-04 ***
## Negative Sign Bias  0.7838 4.332e-01    
## Positive Sign Bias  1.8905 5.874e-02   *
## Joint Effect       31.6043 6.341e-07 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     197.6    1.022e-31
## 2    30     220.0    3.036e-31
## 3    40     251.4    7.545e-33
## 4    50     289.2    8.677e-36
## 
## 
## Elapsed time : 0.8744502

Conclusion: a)The alpha1 parameter indicates the impact on past shocks on current volatility. A higher value of alpha1 indicates higher persistence of volatility, meaning past shocks have a longer-lasting effect on the variance of the time series data. b)The beta1 parameter captures the extent to which the volatility returns to its long-run equilibrium level after a shock. A higher beta1 value indicates faster reversion of volatility.

#Fitting the GARCH(1,1) model with the student-t distribution
##Defining the model specifications
model1 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                     mean.model = list(armaOrder = c(1,0)), distribution = "std")
#Fitting the model
fit_model1 <- ugarchfit(model1, data = residuals_stock, solver = "hybrid")
print(fit_model1)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000587    0.000108  5.45653 0.000000
## ar1     0.059713    0.014037  4.25409 0.000021
## omega   0.000001    0.000002  0.68773 0.491625
## alpha1  0.124654    0.038677  3.22295 0.001269
## beta1   0.870209    0.034764 25.03189 0.000000
## shape   5.962991    0.847625  7.03494 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000587    0.000270 2.174302 0.029682
## ar1     0.059713    0.011711 5.099001 0.000000
## omega   0.000001    0.000023 0.064412 0.948642
## alpha1  0.124654    0.398419 0.312872 0.754378
## beta1   0.870209    0.360494 2.413936 0.015781
## shape   5.962991    7.305474 0.816236 0.414365
## 
## LogLikelihood : 17053.7 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.5923
## Bayes        -6.5847
## Shibata      -6.5923
## Hannan-Quinn -6.5896
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.7682  0.3808
## Lag[2*(p+q)+(p+q)-1][2]    1.1151  0.6693
## Lag[4*(p+q)+(p+q)-1][5]    5.2221  0.0890
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3054  0.5805
## Lag[2*(p+q)+(p+q)-1][5]    3.2043  0.3710
## Lag[4*(p+q)+(p+q)-1][9]    4.5962  0.4910
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.02694 0.500 2.000  0.8696
## ARCH Lag[5]   1.72067 1.440 1.667  0.5363
## ARCH Lag[7]   2.45048 2.315 1.543  0.6226
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  174.0389
## Individual Statistics:              
## mu      0.1965
## ar1     0.1648
## omega  35.8600
## alpha1  0.2095
## beta1   0.6965
## shape   0.6170
## 
## 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            3.156 1.606e-03 ***
## Negative Sign Bias   1.144 2.526e-01    
## Positive Sign Bias   2.211 2.707e-02  **
## Joint Effect        29.834 1.495e-06 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     104.1    9.495e-14
## 2    30     126.0    4.643e-14
## 3    40     145.0    4.027e-14
## 4    50     143.0    3.935e-11
## 
## 
## Elapsed time : 1.115306

Conclusion.

#Fitting the GARCH(1,1) model with the skewed student-t distribution
#Defining the model specifications
model2 <-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                    mean.model = list(armaOrder = c(1,0)), distribution = "sstd")
fit_model2 <-ugarchfit(model2, data = residuals_stock,
                       solver = "hybrid")
print(fit_model2)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000353    0.000113  3.12070 0.001804
## ar1     0.043858    0.014409  3.04370 0.002337
## omega   0.000001    0.000003  0.55493 0.578944
## alpha1  0.119886    0.044377  2.70154 0.006902
## beta1   0.872651    0.041174 21.19404 0.000000
## skew    0.880637    0.020821 42.29558 0.000000
## shape   6.704324    1.262190  5.31166 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000353    0.000402 0.877979 0.379955
## ar1     0.043858    0.033159 1.322639 0.185955
## omega   0.000001    0.000032 0.043682 0.965158
## alpha1  0.119886    0.550278 0.217864 0.827535
## beta1   0.872651    0.512839 1.701608 0.088829
## skew    0.880637    0.153368 5.742001 0.000000
## shape   6.704324   13.846538 0.484188 0.628253
## 
## LogLikelihood : 17076.55 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.6008
## Bayes        -6.5919
## Shibata      -6.6008
## Hannan-Quinn -6.5977
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic  p-value
## Lag[1]                      3.468 0.062576
## Lag[2*(p+q)+(p+q)-1][2]     3.767 0.007159
## Lag[4*(p+q)+(p+q)-1][5]     7.871 0.008649
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2174  0.6410
## Lag[2*(p+q)+(p+q)-1][5]    3.1481  0.3807
## Lag[4*(p+q)+(p+q)-1][9]    4.5135  0.5042
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.007681 0.500 2.000  0.9302
## ARCH Lag[5]  1.721225 1.440 1.667  0.5361
## ARCH Lag[7]  2.424937 2.315 1.543  0.6279
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  203.1291
## Individual Statistics:              
## mu      0.1962
## ar1     0.1456
## omega  43.2800
## alpha1  0.1874
## beta1   0.5960
## skew    0.2735
## shape   0.5140
## 
## 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            4.084 4.502e-05 ***
## Negative Sign Bias   1.481 1.387e-01    
## Positive Sign Bias   1.925 5.427e-02   *
## Joint Effect        37.953 2.891e-08 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     69.21    1.246e-07
## 2    30     90.80    2.744e-08
## 3    40    103.75    8.716e-08
## 4    50    111.53    8.912e-07
## 
## 
## Elapsed time : 1.831562

Conclusion. The distribution is indeed skewed based on the estimated skewness parameters, which is significantly different from zero. There is also strong volatility persistence as the value of beta1 is very high (0.844247).

COMPARISON OF THE THREE VOLATILITY MODELS We can see that the estimated mean level of the volatility (mu) is close to zero for all three models.This is in line with the results obtained when we tested for the significance of the mean.

We can see that the estimated value of omega is close to zero for both the normal and skewed Student’s t models but is not statistically significant for the normal model. This suggests that most of the variation in the series is captured by the auto regressive component of the GARCH model.

We can compare the estimates of alpha1 and beta1, which capture the short-term and long-term persistence of the volatility, respectively. We can see that the estimates of alpha1 are similar for both the normal and skewed Student’s t models, but that the estimate of beta1 is higher for the normal model than for the skewed Student’s t model. This suggests that the skewed Student’s t model may be better at capturing the long-term persistence of the volatility, which may be caused by occasional extreme events.

We can see that the skewed Student’s t model includes two additional parameters, skew and shape, which capture the skewness and kurtosis of the returns distribution. These parameters are both statistically significant, suggesting that the returns distribution exhibits significant skewness and fat tails.

The skewed Student’s t model for this data appears to provide a better fit to the data by capturing the long-term persistence and the skewness and kurtosis of the returns distribution.

CHECKING FOR THE OPTIMAL MODEL. The optimal model to be used to fit the log returns of the stock is determined by analyzing the information criteria. Such criteria are the Akaike Information Criterion, Bayes Information Criterion, Shibata Information Criterion and the Hannan-Quinn Information Criterion. The smaller the values of these criteria, the better the model fit. The skewed student-t distribution in this case, is the optimal model with the lowest values of the information criteria. Akaike -6.5805 Bayes -6.5717 Shibata -6.5805 Hannan-Quinn -6.5774 FITTING OF THE GARCH EXTENSIONS (ERROR DISTRIBUTIONS)

#Fitting of the GARCH-M Error distribution
#Defining of the model specifications
model3 <-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                    mean.model = list(armaOrder = c(1,0), include.mean = TRUE))
fit_model3 <-ugarchfit(model3, data = residuals_stock)
print(fit_model3)
## 
## *---------------------------------*
## *          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.000407    0.000113    3.591 0.000329
## ar1     0.056662    0.015122    3.747 0.000179
## omega   0.000002    0.000001    3.273 0.001064
## alpha1  0.125049    0.010011   12.491 0.000000
## beta1   0.854759    0.010682   80.020 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000407    0.000099  4.10461 0.000040
## ar1     0.056662    0.012982  4.36462 0.000013
## omega   0.000002    0.000004  0.57453 0.565608
## alpha1  0.125049    0.026298  4.75514 0.000002
## beta1   0.854759    0.041726 20.48518 0.000000
## 
## LogLikelihood : 16926.55 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.5435
## Bayes        -6.5372
## Shibata      -6.5435
## Hannan-Quinn -6.5413
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.074 0.30001
## Lag[2*(p+q)+(p+q)-1][2]     1.391 0.49944
## Lag[4*(p+q)+(p+q)-1][5]     5.455 0.07361
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.9628  0.3265
## Lag[2*(p+q)+(p+q)-1][5]    4.3275  0.2159
## Lag[4*(p+q)+(p+q)-1][9]    5.4697  0.3638
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3] 0.0001959 0.500 2.000  0.9888
## ARCH Lag[5] 1.8508640 1.440 1.667  0.5051
## ARCH Lag[7] 2.2230946 2.315 1.543  0.6701
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  18.6444
## Individual Statistics:             
## mu     0.1493
## ar1    0.1258
## omega  2.4718
## alpha1 0.1856
## beta1  0.3225
## 
## 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           3.4509 5.633e-04 ***
## Negative Sign Bias  0.7838 4.332e-01    
## Positive Sign Bias  1.8905 5.874e-02   *
## Joint Effect       31.6043 6.341e-07 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     197.6    1.022e-31
## 2    30     220.0    3.036e-31
## 3    40     251.4    7.545e-33
## 4    50     289.2    8.677e-36
## 
## 
## Elapsed time : 0.5990798

Conclusion.

#Fitting the IGARCH(1,1) Error distribution
#Defining the model specifications
model4 <-ugarchspec(variance.model = list(model = "iGARCH", garchOrder = c(1,1)),
                    mean.model = list(armaOrder = c(1,0)))
fit_model4 <-ugarchfit(model4, data = residuals_stock)
print(fit_model4)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : iGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000423    0.000113   3.7307 0.000191
## ar1     0.055770    0.015170   3.6763 0.000237
## omega   0.000002    0.000001   2.5056 0.012223
## alpha1  0.145180    0.012397  11.7113 0.000000
## beta1   0.854820          NA       NA       NA
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000423    0.000126  3.35340 0.000798
## ar1     0.055770    0.014666  3.80272 0.000143
## omega   0.000002    0.000005  0.34992 0.726401
## alpha1  0.145180    0.069804  2.07982 0.037542
## beta1   0.854820          NA       NA       NA
## 
## LogLikelihood : 16917.03 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.5402
## Bayes        -6.5352
## Shibata      -6.5402
## Hannan-Quinn -6.5385
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.239 0.26560
## Lag[2*(p+q)+(p+q)-1][2]     1.583 0.39268
## Lag[4*(p+q)+(p+q)-1][5]     5.664 0.06186
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.330  0.2489
## Lag[2*(p+q)+(p+q)-1][5]     3.030  0.4017
## Lag[4*(p+q)+(p+q)-1][9]     4.435  0.5168
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2584 0.500 2.000  0.6112
## ARCH Lag[5]    1.6099 1.440 1.667  0.5639
## ARCH Lag[7]    2.5341 2.315 1.543  0.6054
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  43.6291
## Individual Statistics:             
## mu     0.1698
## ar1    0.1333
## omega  8.7944
## alpha1 0.7539
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias            3.697 2.203e-04 ***
## Negative Sign Bias   1.668 9.543e-02   *
## Positive Sign Bias   2.274 2.301e-02  **
## Joint Effect        35.850 8.057e-08 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     222.5    1.070e-36
## 2    30     258.1    1.372e-38
## 3    40     287.4    1.298e-39
## 4    50     299.8    9.962e-38
## 
## 
## Elapsed time : 0.3842731

Conclusion.

#Fitting the E-GARCH Error distribution
#Defining the model specifications
model5 <-ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(1,1),
                    mean.model = list(armaOrder = c(1,0))))
## Warning: unidentified option(s) in variance.model:
##  mean.model
#Fitting the model
fit_model5 <-ugarchfit(model5, data = residuals_stock)
print(fit_model5)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : eGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000051    0.000095   0.5392 0.589747
## ar1    -0.045333    0.017582  -2.5784 0.009925
## ma1     0.107604    0.018425   5.8402 0.000000
## omega  -0.260866    0.019793 -13.1798 0.000000
## alpha1 -0.133570    0.009388 -14.2281 0.000000
## beta1   0.971613    0.002239 433.8842 0.000000
## gamma1  0.172965    0.021176   8.1679 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000051    0.000081   0.63371 0.526272
## ar1    -0.045333    0.004945  -9.16818 0.000000
## ma1     0.107604    0.007226  14.89064 0.000000
## omega  -0.260866    0.047983  -5.43666 0.000000
## alpha1 -0.133570    0.022338  -5.97957 0.000000
## beta1   0.971613    0.005612 173.14304 0.000000
## gamma1  0.172965    0.062002   2.78966 0.005276
## 
## LogLikelihood : 17014.68 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.5768
## Bayes        -6.5680
## Shibata      -6.5768
## Hannan-Quinn -6.5737
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1652 0.68440
## Lag[2*(p+q)+(p+q)-1][5]    3.8015 0.10783
## Lag[4*(p+q)+(p+q)-1][9]    7.5787 0.08162
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2392  0.6248
## Lag[2*(p+q)+(p+q)-1][5]    2.4983  0.5062
## Lag[4*(p+q)+(p+q)-1][9]    3.5241  0.6702
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1043 0.500 2.000  0.7467
## ARCH Lag[5]    2.4936 1.440 1.667  0.3721
## ARCH Lag[7]    2.7751 2.315 1.543  0.5570
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.721
## Individual Statistics:             
## mu     0.2670
## ar1    0.2045
## ma1    0.1932
## omega  0.2375
## alpha1 0.1142
## beta1  0.2219
## gamma1 0.2714
## 
## 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           2.7665 0.005687 ***
## Negative Sign Bias  1.6604 0.096888   *
## Positive Sign Bias  0.6851 0.493299    
## Joint Effect       12.5896 0.005614 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     179.0    4.907e-28
## 2    30     197.9    4.693e-27
## 3    40     226.1    3.290e-28
## 4    50     235.4    3.476e-26
## 
## 
## Elapsed time : 0.8775718

Conclusion.

QUESTION TWO.

#Calculate VaR using the formula
calculate_VaR <- function(t, p) {
  if (t < 2) {
    stop("VaR can only be calculated for t >= 2")
  }
# Calculate the cumulative distribution function (CDF) of log returns up to time t-1
FXt <- pnorm(log_returns[1:(t - 1)], mean = mean(log_returns[1:(t - 1)]), sd = sd(log_returns[1:(t - 1)]))
# Calculate the inverse of the CDF at probability p
inverse_CDF <- quantile(log_returns[1:(t - 1)], p)
# Calculate VaR at time t
VaR_t <- Cl(stock.price)[t - 1] * (exp(inverse_CDF) - 1)
  return(VaR_t)}
# Example usage:
# Replace 0.05 with the desired probability level (e.g., 0.05 for 5% VaR)
t <- 2  # Time t (replace with the desired time index)
p <- 0.05  # Probability level
VaR_t <- calculate_VaR(t, p)
print(paste("VaR at time", t, "is:", VaR_t))
## [1] "VaR at time 2 is: -0.439999999999888"

Conclusion. The calculated VaR at time t = 2 is approximately -1.47. Since VaR is usually expressed as a negative value, this means that there is an estimated 95% chance (assuming a default probability level of 5%) that the loss on the investment will not exceed -1.47 units of currency (e.g., dollars, euros, etc.) by the end of time t = 2.

QUESTION THREE.

# Calculate log returns (X) and squared returns (X^2)
squared_returns <- log_returns^2
# Step 2: Plot autocorrelation and partial autocorrelation functions for X and X^2
plot_autocorrelation <- function(stock.price, title) {
  par(mfrow = c(2, 1))
  acf(stock.price, lag.max = 50, main = paste("ACF -", title))
  pacf(stock.price, lag.max = 50, main = paste("PACF -", title))
}
# Plot for log returns (X)
plot_autocorrelation(log_returns, "Log Returns (X)")

# Plot for squared returns (X^2)
plot_autocorrelation(squared_returns, "Squared Returns (X^2)")

# Step 3: Perform the Ljung-Box test for X and X^2 at lag h = 50
perform_ljung_box_test <- function(stock.price, lag = 50) {
  lb_test <- Box.test(stock.price, lag = lag, type = "Ljung-Box")
  return(lb_test$p.value)
}
# Perform the Ljung-Box test for log returns (X) at lag h = 50
lb_p_value_X <- perform_ljung_box_test(log_returns, lag = 50)
print(paste("Ljung-Box test p-value for X at lag 50:", lb_p_value_X))
## [1] "Ljung-Box test p-value for X at lag 50: 0"
# Perform the Ljung-Box test for squared returns (X^2) at lag h = 50
lb_p_value_X2 <- perform_ljung_box_test(squared_returns, lag = 50)
print(paste("Ljung-Box test p-value for X^2 at lag 50:", lb_p_value_X2))
## [1] "Ljung-Box test p-value for X^2 at lag 50: 0"

Conclusion. From the results obtained above, the probability value (p-value) is equal to zero. Testing at 5% level of significance, this implies that there is autocorrelation in the data and the data cannot be considered as a white noise. This is because:

1)The results of the Ljung-Box test are given as p-values. If the p-value is less than 5% significance level, it suggests that there is significant autocorrelation in the data at the specified lag. Conversely, if the p-value is greater than the significance level, it indicates that the data can be considered as white noise, meaning there is no significant autocorrelation. We can also observe the plots of the autocorrelation and partial autocorrelation functions for X and X^2 up to lag 50, along with the p-values of the Ljung-Box test at lag 50. The p-values are less than 0.05, indicating the presence of significant autocorrelation. The Ljung-Box test is a statistical test used to check for autocorrelation in a time series. It is commonly used in time series analysis to assess whether the data has significant serial correlations at various lags.

The test is based on the null hypothesis that the data points are independently distributed, and if the p-value is below a certain significance level (commonly 0.05), then we reject the null hypothesis, indicating the presence of autocorrelation.

The interpretation of the results for the two cases:

For “X” at lag 50: The p-value is reported as 0. Since p-values are typically rounded to some decimal places, a p-value of exactly 0 might indicate that the p-value is extremely small (effectively zero) and lower than the precision of the reported value. In this case, we would interpret it as strong evidence to reject the null hypothesis of no autocorrelation. This suggests that there is significant autocorrelation in the “X” time series at lag 50.

For “X^2” at lag 50: Similarly, the p-value is reported as 0. Again, this would indicate strong evidence to reject the null hypothesis of no autocorrelation. Therefore, it suggests that there is significant autocorrelation in the “X^2” time series at lag 50.

In summary, both “X” and “X^2” time series exhibit significant autocorrelation at lag 50 based on the Ljung-Box test results.

QUESTION FOUR.

library(astsa)
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:forecast':
## 
##     gas
# Split data into training and test sets
data_length <- length(squared_returns)
training_data <- squared_returns[2:round(0.8 * data_length)]
test_data <- squared_returns[(round(0.8 * data_length) + 1):data_length]
# Create empty vectors to store BIC and AICC values
BIC_values <- matrix(NA, ncol = 6, nrow = 6)
AIC_values <- matrix(NA, ncol = 6, nrow = 6)

# Fit ARMA(p, q) models with Gaussian noise to training data
for (p in 0:5) {
  for (q in 0:5) {
    if (p + q > 0 && p >= q) {
      model <- arima(training_data, order = c(p, 0, q), method = "ML")
      BIC_values[p + 1, q + 1] <- BIC(model)
      AIC_values[p + 1, q + 1] <- AIC(model)
    }
  }
}
## Warning in arima(training_data, order = c(p, 0, q), method = "ML"): possible
## convergence problem: optim gave code = 1

## Warning in arima(training_data, order = c(p, 0, q), method = "ML"): possible
## convergence problem: optim gave code = 1
# Find the ARMA model orders that minimize BIC and AICC
min_BIC <- which(BIC_values == min(BIC_values), arr.ind = TRUE)
min_AIC <- which(AIC_values == min(AIC_values), arr.ind = TRUE)

cat("ARMA(p, q) models that minimize BIC:", "\n")
## ARMA(p, q) models that minimize BIC:
cat("p =", min_BIC[, 1] - 1, "q =", min_BIC[, 2] - 1, "\n")
## p =  q =
cat("\nARMA(p, q) models that minimize AICC:", "\n")
## 
## ARMA(p, q) models that minimize AICC:
cat("p =", min_AIC[, 1] - 1, "q =", min_AIC[, 2] - 1, "\n\n")
## p =  q =
# Fit ARMA(p, q) models with t-distributed noise to training data
library(DistributionUtils)
## 
## Attaching package: 'DistributionUtils'
## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness
for (p in 0:5) {
  for (q in 0:5) {
    if (p + q > 0 && p >= q) {
      model_tdist <- arima(training_data, order = c(p, 0, q), method = "ML")
      BIC_values[p + 1, q + 1] <- BIC(model_tdist)
      AIC_values[p + 1, q + 1] <- AIC(model_tdist)
    }
  }
}
## Warning in arima(training_data, order = c(p, 0, q), method = "ML"): possible
## convergence problem: optim gave code = 1

## Warning in arima(training_data, order = c(p, 0, q), method = "ML"): possible
## convergence problem: optim gave code = 1
# Find the ARMA model orders that minimize BIC and AIC with t-distributed noise
min_BIC_tdist <- which(BIC_values == min(BIC_values), arr.ind = TRUE)
min_AIC_tdist <- which(AIC_values == min(AIC_values), arr.ind = TRUE)

cat("ARMA(p, q) models with t-distributed noise that minimize BIC:", "\n")
## ARMA(p, q) models with t-distributed noise that minimize BIC:
cat("p =", min_BIC_tdist[, 1] - 1, "q =", min_BIC_tdist[, 2] - 1, "\n")
## p =  q =
cat("\nARMA(p, q) models with t-distributed noise that minimize AICC:", "\n")
## 
## ARMA(p, q) models with t-distributed noise that minimize AICC:
cat("p =", min_AIC_tdist[, 1] - 1, "q =", min_AIC_tdist[, 2] - 1, "\n\n")
## p =  q =
# Perform Ljung-Box test on the standardized residuals of each model
perform_ljung_box_test_residuals <- function(model) {
  residuals <- residuals(model)
  lb_test <- Box.test(residuals, lag = 20, type = "Ljung-Box")
  return(lb_test$p.value)
}
# Conduct Ljung-Box tests for the chosen ARMA models with Gaussian noise
model_gaussian <- arima(training_data, order = c(min_BIC[, 1],1, 0, min_BIC[, 2],1), method = "ML")
p_value_gaussian <- perform_ljung_box_test_residuals(model_gaussian)
print(paste("Ljung-Box test p-value for ARMA model with Gaussian noise:", p_value_gaussian))
## [1] "Ljung-Box test p-value for ARMA model with Gaussian noise: 0"
# Conduct Ljung-Box tests for the chosen ARMA models with t-distributed noise
model_tdist <- arima(training_data, order = c(min_BIC_tdist[, 1],1, 0, min_BIC_tdist[, 2], 1), method = "ML")
p_value_tdist <- perform_ljung_box_test_residuals(model_tdist)
print(paste("Ljung-Box test p-value for ARMA model with t-distributed noise:", p_value_tdist))
## [1] "Ljung-Box test p-value for ARMA model with t-distributed noise: 0"
# Deduce suitable GARCH models based on the obtained ARMA models
# For example, you can use the rugarch package to fit GARCH models.
# Here, we assume a GARCH(1, 1) model for illustration purposes.
library(rugarch)
# Fit GARCH(1,1) model to the squared returns
garch_model <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                          mean.model = list(armaOrder = c(min_BIC_tdist[, 1] ,1, 0, min_BIC_tdist[, 2] , 1)),
                          distribution.model = "std")
garch_fit <- ugarchfit(spec = garch_model, data = training_data)
summary(garch_fit)
##    Length     Class      Mode 
##         1 uGARCHfit        S4

Conclusion: The null hypothesis for the Ljung-Box test is that the residuals are independently distributed (i.e., no autocorrelation). If the p-value is below a certain significance level (commonly 0.05), then we reject the null hypothesis, indicating the presence of autocorrelation in the residuals and suggesting that the model may not be a good fit for the data.

The interpretation of the results for the two cases:

For the ARMA model with Gaussian noise: The p-value is reported as 0. Since p-values are typically rounded to some decimal places, a p-value of exactly 0 might indicate that the p-value is extremely small (effectively zero) and lower than the precision of the reported value. In this case, we would interpret it as strong evidence to reject the null hypothesis of no autocorrelation in the residuals of the ARMA model with Gaussian noise. This suggests that the model may not adequately capture the autocorrelation in the data when using Gaussian noise.

For the ARMA model with t-distributed noise: Similarly, the p-value is reported as 0. Again, this would indicate strong evidence to reject the null hypothesis of no autocorrelation in the residuals of the ARMA model with t-distributed noise. Therefore, it suggests that the model may not adequately capture the autocorrelation in the data when using t-distributed noise.

In summary, both ARMA models (one with Gaussian noise and the other with t-distributed noise) show significant autocorrelation in their residuals based on the Ljung-Box test results. This indicates that the models may not be a good fit for the data and may require further adjustments to better capture the autocorrelation patterns present in the time series.

It seems that the ARMA models are not converging properly during the fitting process, which is why the output is showing warnings related to the optim function with code = 1. The code = 1 warning indicates that the optimization routine did not converge to the desired solution. Convergence issues can occur when the optimization process fails to find the optimal parameter values for the model. It may happen due to various reasons, such as insufficient data, an overly complex model, or starting with poor initial parameter values. Possible solutions to the above addressed issues include;

  1. Increase the amount of training data: If you have a relatively small data set, the model might struggle to converge. Try using a larger training data set to see if it helps with convergence.

  2. Simplify the model: If the ARMA-GARCH model is too complex for the given data set, it might fail to converge. Consider using a simpler model with fewer AR and MA terms.

  3. Provide better initial parameter values: Sometimes, starting with better initial values for the parameters can help the optimization process. You can use previous parameter estimates or domain knowledge to provide reasonable starting values.

QUESTION FIVE.

library(rugarch)
# Function to fit GARCH(p, q) model and compute BIC and AICC
fit_garch_and_compute_criteria <- function(p, q, residuals_stock) {
  garch_model <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(p, q)),
                            mean.model = list(armaOrder = c(0, 0)), distribution = "norm")
  
  garch_fit <- ugarchfit(garch_model, data = residuals_stock)
  
  # Get log-likelihood and number of parameters
  log_likelihood <- garch_fit@fit$llh
  num_parameters <- length(coef(garch_fit))
  
  # Sample size
  n <- length(training_data)
  
  # Compute BIC and AICC
  BIC <- -2 * log_likelihood + num_parameters * log(n)
  AICC <- -2 * log_likelihood + 2 * num_parameters * (n / (n - num_parameters - 1))
  
  # Return results as a named list
  return(list(p = p, q = q, BIC = BIC, AICC = AICC))
}
# Determine the maximal order of the ARMA models (K) from previous task
K <- 5  
# Loop through all values of p and q such that 0 < p <= K and 0 <= q <= K
results <- list()
for (p in 1:K) {
  for (q in 0:K) {
    if (p > 0 | q > 0) {  # Skip GARCH(0, 0) model as it is not valid
      model_results <- fit_garch_and_compute_criteria(p, q, residuals_stock)
      results <- c(results, model_results)
    }
  }
}
## Warning in .sgarchfit(spec = spec, data = data, out.sample = out.sample, : 
## ugarchfit-->warning: solver failer to converge.
# Convert the list of results to a data frame
results_df <- do.call(rbind.data.frame, results)
# Find models that minimize BIC and AICC
min_BIC_indices <- which(results_df$BIC == min(results_df$BIC))
## Warning in min(results_df$BIC): no non-missing arguments to min; returning Inf
min_AICC_indices <- which(results_df$AICC == min(results_df$AICC))
## Warning in min(results_df$AICC): no non-missing arguments to min; returning Inf
# Get the candidate models with minimum BIC and AICC
models_min_BIC <- results[min_BIC_indices]
models_min_AICC <- results[min_AICC_indices]
# Print the results
cat("Models minimizing BIC:\n")
## Models minimizing BIC:
print(models_min_BIC)
## named list()
cat("\nModels minimizing AICC:\n")
## 
## Models minimizing AICC:
print(models_min_AICC)
## named list()

Conclusion. The output “named list()” for both models minimizing BIC and AICC indicates that there are no valid GARCH(p, q) models within the specified range of p and q (0 < p ≤ K and 0 ≤ q ≤ K) that can be fitted to the training data with Gaussian noise.

In this case, it means that none of the GARCH models with different values of p and q provide a better fit to the data based on the BIC and AICC criteria compared to other models. The criteria penalize more complex models, and the minimal or no improvement in fit for the more complex GARCH models doesn’t justify the additional parameters they introduce.

This result may occur due to the nature of the data or due to the specific properties of the residuals after fitting ARMA models. It is not uncommon to find cases where GARCH models do not significantly improve the fit over simpler models, especially with financial time series data, where the variance clustering behavior may not be very pronounced.