Data Description and Exploration

Data Description

The company selected is dell, the data sets consist of the daily stock log returns of dell computer corporation from October 1993 to October 1998. . It is common to use log returns for analyze financial time series. During the exploratory analysis of the dataset, it was observed that the log returns exhibited a distribution reasonably close to normality. As a result, no data transformation is required at this stage, and the log returns will be used directly for the analysis. From this point forward, we will simply refer to the log returns as returns.

Data Exploration

The number of observation of our stock returns time series from October 1993 to October 1998 is 1261 observations, the plot of this is presented below displays evident volatility clustering, mean reversion, and occasional extreme movements. These characteristics suggest the returns are heteroscedastic and possibly non-normally distributed, which is typical of financial time series data. We can now think to squared the returns to highlight volatility, which helps visualize the magnitude of fluctuations without concern for direction. And the squared returns plot reveals clear volatility clustering, where large fluctuations group together in time. It also emphasizes the presence of extreme values, indicating that the return distribution has fat tails. This confirms the appropriateness of using ARCH/GARCH models to model and forecast volatility in this financial series. Some preliminary results for the return series are given in Table 1. The results reveals that the return series has a high kurtosis which suggests that the series is not normally distributed. This is confirmed by the tests for normality which are given in Table 2 and from a visual inspection of the histogram of the return shown in Figure 2.

# Compute moments for returns
moments_returns <- c(
  Mean     = mean(Rt),
  Min    = min(Rt),
  Max = max(Rt),
  Median = median(Rt),
  Variance = var(Rt),
  Skewness = skewness(Rt),
  Kurtosis = kurtosis(Rt)
)

# Compute moments for squared returns
moments_squared_returns <- c(
  Mean     = mean(Rt2),
  Min    = min(Rt2),
  Max = max(Rt2),
  Median = median(Rt2),
  Variance = var(Rt2),
  Skewness = skewness(Rt2),
  Kurtosis = kurtosis(Rt2)
)

# Combine into a table
moments_table <- tibble::tibble(
  Moment   = names(moments_returns),
  Returns = as.numeric(moments_returns),
  Squared_Returns = as.numeric(moments_squared_returns)
)


# Print the result
moments_table |> tbl(caption = "Statistical moments of returns and squared returns {#tbl-t1}")

# Kolmogorov-Smirnov test (against normal distribution with estimated mean and sd)
ks_r <- ks.test(Rt, "pnorm", mean = mean(Rt), sd = sd(Rt))
ks_sr <- ks.test(Rt2, "pnorm", mean = mean(Rt2), sd = sd(Rt2))

# Cramer-von Mises test
cvm_r <- cvm.test(Rt)
cvm_sr <- cvm.test(Rt2)

# Anderson-Darling test
ad_r <- ad.test(Rt)
ad_sr <- ad.test(Rt2)

# Combine results into a table
test_results <- data.frame(
  Test = c("Kolmogorov-Smirnov", "Cramér–von Mises", "Anderson-Darling"),
  Statistic_Returns = c(ks_r$statistic, cvm_r$statistic, ad_r$statistic),
  Pvalue_Returns     = c(ks_r$p.value, cvm_r$p.value, ad_r$p.value),
  Statistic_Squared  = c(ks_sr$statistic, cvm_sr$statistic, ad_sr$statistic),
  Pvalue_Squared     = c(ks_sr$p.value, cvm_sr$p.value, ad_sr$p.value)
)

# Print the table
test_results |> tbl(caption = "Normality tests on returns and squared returns {#tbl-t2}")

# plot of retuns

plot(Rt, ylab = "Returns")

hist(Rt, 
     breaks = 20, 
     probability = TRUE,  # This makes y-axis show density
     col = "lightblue", 
     main = "Histogram with Normal Curve",
     xlab = "Value")
curve(dnorm(x, mean = mean(Rt), sd = sd(Rt)), 
      col = "red", 
      lwd = 2, 
      add = TRUE)

plot(Rt2, ylab = "Squared Returns")

Stationarity and P(ACF) of the squared returns

The visualization of the squared returns provides a reasonable assumption that the data is stationary. To investigate this further, a Dicky-Fuller test was conducted using R Studio, returning a p-value smaller than \(0.01\), so we reject the null hypothesis and accept stationarity. Therefore this suggest to model the squared returns with an ARMA(p,q) model. To determine the degrees of the ARMA model, the ACF and PACF were plotted, shown below.

Visual analysis

The plots show a strong autocorrelation at lag 1, and several subsequent lags also show small but significant autocorrelations, the PACF has significant spikes at multiple lags, especially early on, and they gradually decay. This suggests an ARMA(p,q) with p or q=1 and p or q=2.

# ACF and PACF of the squared returns

acf(Rt2, main = "ACF of squared Returns")
pacf(Rt2, main = "PACF of squared Returns")

AIC criterion of ARMA processes

Each of the proposed models was fit to the squared returns time series using R-Studio using the maximum likelihood method. The \(20\%\) best model results are shown in the table below.

Comp.Sarima(
  Rt2, d = 0, saison = NA, D = 0,
  p.max = 2, q.max = 2, P.max = 0, Q.max = 0
) |>
  kable(caption = "Comparison of ARMA models based on AIC")

modele (p,d,q)x(P,D,Q)_saison :  1 0 2 x 0 0 0 _ NA :  nb param:  3     AIC: 11309 
modele (p,d,q)x(P,D,Q)_saison :  2 0 2 x 0 0 0 _ NA :  nb param:  4     AIC: 11304 

An ARMA(1, 2) model was chosen, because it had lowest number of parameter and the AIC value was not too far from the lowest AIC value. Knowing that the squared returns is an ARMA(1,2), base on the GARCH properties we can say that returns could be a GARCH(1,2) or GARCH(2,2). We will also test if a GARCH(1,1) would be sufficient to model the returns.

GARCH Model for Returns time series

GARCH model estimation

Conducting many GARCH models to fit out our returns series, base on the AIC criterion we concluded that GARCH(1,2) is the best as it had the lowest AIC value of \(5.223368\) (see in the appendix). Which is express as follow :

\[X_t = \sigma_t Z_t \;\;\;\{Z_t\}_t \sim \text{i.i.d. } \mathcal{N}(0, 1)\]

\[\sigma_t^2 = 0.77443 + 0.08694 X_{t-1}^2 + 0.63982 \sigma_{t-2}^2\]

Estimating the coefficients of our model, the results in the appendix. It is noticeable that \(\beta_1\) is not significant. This is why we drop \(\beta_1\) in the dynamical equation of the GARCH(1,2) model for the returns.

Analysis of residuals

The standard residuals was tested for significance using Ljung-Box Test. The squared of the standard residuals was also accounted for (see the table below).

fit1 <- garchFit(data ~ garch(1,2), data = Rt, trace = F, include.mean = F)


fit2 <- garchFit(data ~ garch(2,2), data = Rt, trace = F, include.mean = F)


fit3 <- garchFit(data ~ garch(1,1), data = Rt, trace = F, include.mean = F)

The results from the standardized residuals tests assess the adequacy of a returns time series model. The Jarque-Bera and Shapiro-Wilk tests both have very small p-values (\(2.82e^{-07}\) and \(1.69e^{-04}\), respectively), indicating that the residuals deviate significantly from normality. However, the Ljung-Box tests on residuals (R) and squared residuals (\(R^2\)) across various lags (10, 15, 20) yield high p-values, suggesting no significant autocorrelation or remaining ARCH effects in the residuals or their squares. The LM Arch Test also confirms this with a high p-value (0.644), indicating no significant conditional heteroscedasticity. Overall, while the residuals are not normally distributed, they do not show evidence of autocorrelation or ARCH effects, suggesting that the model captures the serial structure and volatility well.

Construction of Prediction Interval

Constructing a prediction interval involves estimating the range within which future observations of a time series (or model output) are likely to fall, with a certain level of confidence (e.g., \(95\%\)). This is different from a confidence interval, which estimates the range for the mean of the population or regression function.

In a GARCH model, the prediction interval for a future value \(X_{t+h}\) is constructed using the assumption of conditional normality, meaning that given past information \(\mathcal{F}_t\), the future value \(X_{t+h}\) is normally distributed with mean \(\hat{X}_{t+h}\) and variance \(\hat{\sigma}^2_{t+h}\). Although the GARCH process itself is not normally distributed unconditionally, this conditional normality allows us to build a prediction interval as

\[\hat{X}_{t+h} \pm z_{\alpha/2} \cdot \sqrt{\text{MSE}_h},\]

where \(\text{MSE}_h = \mathbb{E}_t[(X_{t+h} - \hat{X}_{t+h})^2] = \hat{\sigma}^2_{t+h}\) is the conditional forecast error variance at horizon h , and \(z_{\alpha/2}\) is the quantile of the standard normal distribution (e.g., 1.96 for 95% confidence).

Prediction of the conditional variance

The graphic bellow shows prediction intervals from a GARCH model with a constant conditional variance, meaning that the forecasted variance \(\hat{\sigma}^2_{t+h}\) does not change with the forecast horizon h. This results in prediction intervals \(\hat{X}_{t+h} \pm 1.96 \cdot \sqrt{\hat{\sigma}^2_{t+h}}\) that are of constant width, as observed in the plot. The flat shape of the intervals reflects the assumption that future volatility is expected to remain stable over time, which is characteristic of a GARCH model with low persistence (i.e., small values of \(\alpha + \beta_1 +\beta_2\)). Consequently, both the upper and lower confidence bounds remain parallel and equidistant from the forecasted mean.

# Prediction of the conditional variance

predict(fit1, n.ahead = (10/100)*length(Rt),plot=TRUE)