Financial market is dynamic and full of uncertainty, understanding and predicting financial volatility thus becomes a topic that catches a lot of investors’ interest. NASDAQ is a famous stock market index that comprises many companies in the field of information technology.[1] It is a good indicator of the dynamics of financial market volatility. Given its importance, modeling it and making predictions based on the model can be beneficial for understanding the changes in the whole financial market and making wise decisions for the future.
In this project, we intend to figure out the best method to model the volatility of NASDAQ log returns in 5 years. We will fit different models, including ARMA, GARCH, and POMP, to compare their results and identify the most appropriate one.
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## -- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
## v dplyr 1.1.4 v readr 2.1.5
## v forcats 1.0.0 v stringr 1.5.1
## v ggplot2 3.5.1 v tibble 3.2.1
## v lubridate 1.9.3 v tidyr 1.3.1
## v purrr 1.0.2
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
##
## Attaching package: 'tsibble'
##
##
## The following object is masked from 'package:lubridate':
##
## interval
##
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, uniontime_series_data <- read.csv("NDAQ.csv", header = TRUE, stringsAsFactors = FALSE)
time_series_data <- time_series_data %>%
mutate(Log_Return = c(NA, diff(log(Close))))
ts_data <- ts(time_series_data$Log_Return, start = c(2019, 4), frequency = 1)# Plot the time series
plot(ts_data, main = "Time Series", xlab = "Time", ylab = "Value")
The following plot is created for visualizing the distribution of the Log_Return variable, in the time series.
# Histogram of the time series values
hist(ts_data, main = "Histogram of Time Series", xlab = "Value", breaks = 30)
# Boxplot to check for outliers
boxplot(ts_data, main = "Boxplot of Time Series")
The ACF measures the correlation between a time series and its lagged values at different time lags. It helps in identifying the presence of seasonality and determining the order of autoregressive (AR) and moving average (MA) terms in time series models. Peaks or significant spikes in the ACF plot indicate potential seasonal patterns or lagged relationships in the data.
On the other hand, the PACF measures the correlation between a time series and its lagged values, while adjusting for the intermediate lags. It helps in identifying the direct relationship between observations at different time lags, thus assisting in determining the order of the AR terms in time series models. Significant spikes in the PACF plot indicate the number of AR terms needed to adequately model the data.
# Autocorrelation function (ACF) and Partial Autocorrelation function (PACF)
ts_data <- na.omit(ts_data)
acf(ts_data, main = "ACF of Time Series")
pacf(ts_data, main = "PACF of Time Series")
The number of significant spikes in the ACF plot is 1, hence, we can assume that the AR term has value 1. Likewise, the number of significant spikes in the PACF plot is 4. Hence, it can be inferred that the MA term is 4.
In the below provided code, adf.test(ts_data, alternative = "stationary") conducts the ADF test on the time series data stored in ts_data. The argument alternative = “stationary” specifies that the null hypothesis of the test is that the series is non-stationary and the alternative hypothesis is that the series is stationary.
# Augmented Dickey-Fuller Test for stationarity
adf_result <- adf.test(ts_data, alternative = "stationary")
print(adf_result)##
## Augmented Dickey-Fuller Test
##
## data: ts_data
## Dickey-Fuller = -10.246, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationaryThe obtained p-value of 0.01 is smaller than the printed p-value, suggesting strong evidence against the null hypothesis. Therefore, we reject the null hypothesis in favor of the alternative, indicating that the time series is stationary.
# KPSS test for stationarity
kpss_result <- kpss.test(ts_data)
print(kpss_result)##
## KPSS Test for Level Stationarity
##
## data: ts_data
## KPSS Level = 0.134, Truncation lag parameter = 7, p-value = 0.1Since the p-value is greater than the significance level, it suggests that there is not enough evidence to reject the null hypothesis. In other words, the data is stationary.
From both the above test results, it can be inferred that the given time series data is stationary. Hence, the data need not be differenced in order to achieve stationarity.
We first consider constructing an ARMA model on the log returns.
In order to choose the optimal values of \(p\) and \(q\) in the ARMA(p,q) model, we use AIC as a selection criterion.
## Loading required package: knitr| MA0 | MA1 | MA2 | MA3 | MA4 | MA5 | |
|---|---|---|---|---|---|---|
| AR0 | -6576.49 | -6581.54 | -6581.60 | -6592.26 | -6592.66 | -6592.84 |
| AR1 | -6582.13 | -6583.29 | -6582.59 | -6591.62 | -6610.24 | -6608.61 |
| AR2 | -6582.65 | -6582.15 | -6585.67 | -6590.85 | -6608.99 | -6607.72 |
| AR3 | -6588.47 | -6590.55 | -6598.05 | -6598.12 | -6616.52 | -6607.83 |
| AR4 | -6595.49 | -6616.24 | -6614.24 | -6617.54 | -6629.82 | -6627.85 |
| AR5 | -6595.46 | -6614.24 | -6612.57 | -6614.96 | -6627.84 | -6625.83 |
\(p\) and \(q\) values are searched in \(\{0,1,2,3,4,5\} \times \{0,1,2,3,4,5\}\) and AIC values of the corresponding models are calculated. As suggested in the AIC table above, ARMA(4,4) model has the lowest AIC value, and thus should be chosen.
nasdaq_arma44 <- arima(ts_data, order = c(4,0,4))
nasdaq_arma44##
## Call:
## arima(x = ts_data, order = c(4, 0, 4))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## -0.5070 0.4405 -0.4425 -0.7861 0.4533 -0.4350 0.4119 0.6480
## s.e. 0.0617 0.0777 0.0886 0.0659 0.0752 0.0923 0.1036 0.0788
## intercept
## 6e-04
## s.e. 4e-04
##
## sigma^2 estimated as 0.000295: log likelihood = 3324.91, aic = -6629.82This ARMA(4,4) model has a log likelihood of 3324.91.
In the ACF plot, almost all auto-correlation coefficients fall within the confidence boundary. The auto-correlation is only significant at lag 18 and lag 28. The QQ-plot shows that the residuals have significantly heavier tails on both sides than normal distribution. This may suggest that ARMA model is not very suitable for fitting the log returns. Therefore, we should try other models which can give better fit to the data.
GARCH model is often employed in the field of finance since it is good at capturing the volatility in financial time series data.[1] Specifically, we say a process \(X = \{X_n\}\) is GARCH(p,q) if \[ X_n = \mu_n + \sigma_n \epsilon_n \] where \(\{\epsilon_n\}\) is an iid white noise process with mean 0 and variance of 1, and the model for \(\sigma_n\) is \[ \sigma_n^2 = \alpha_0 + \sum_{i=1}^p \beta_i \sigma_{n-i}^2 + \sum_{j=1}^q \alpha_j (X_{n-j} - \mu_{n-j})^2 \]
Firstly, we assume that the white noise process \(\{\epsilon_n\}\) in the GARCH model follows normal distribution. Based on this hypothesis, we conduct model selection and model diagnostics to judge whether this assumption is appropriate.
In order to obtain the optimal value of \(p\) and \(q\) in a GARCH(p,q) model under normal white noise, we use AIC value as the criterion here. The model with the smallest AIC value will be considered as the best.
## 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| q=1 | q=2 | q=3 | q=4 | q=5 | |
|---|---|---|---|---|---|
| p=1 | -5.51138 | -5.50953 | -5.50950 | -5.50738 | -5.50533 |
| p=2 | -5.50955 | -5.51054 | -5.51090 | -5.50878 | -5.50656 |
| p=3 | -5.50747 | -5.50588 | -5.50931 | -5.50719 | -5.50497 |
| p=4 | -5.51719 | -5.51560 | -5.51407 | -5.51311 | -5.51116 |
| p=5 | -5.51541 | -5.51430 | -5.51138 | -5.51162 | -5.51016 |
As indicated in the table above, GARCH(4,1) model has the lowest AIC value, and could possibly be the best model.
likelihood(nasdaq_garch41_normal)## [1] 3476.553log(likelihood(nasdaq_garch41_normal))## [1] 8.153797This model has a likelihood of 3476.553 and a log likelihood of 8.1538.
We further generate a QQ-plot and an ACF plot of the GARCH(4,1) model residuals to observe the distribution and to judge whether the auto-correlations are weak.
Disappointingly, auto-correlations at lag 3, 13, 24, 30 seem to be significant. In addition, the distribution of the residuals still have heavier tails on both sides than the normal distribution. This motivates us to go beyond the normally distributed white noise assumption and study some other distributions that have heavier tails, such as t distribution.
In order to better fit the heavy tails of the log return data, we now assume that the white noise process \(\{\epsilon_n\}\) in the GARCH model follows iid t distribution.[2]
Based on AIC criterion again, we conduct a selection of \(p\) and \(q\) value in GARCH(p,q).
| q=1 | q=2 | q=3 | q=4 | q=5 | |
|---|---|---|---|---|---|
| p=1 | -5.61777 | -5.61657 | -5.61595 | -5.61388 | -5.61188 |
| p=2 | -5.61595 | -5.61498 | -5.61469 | -5.61264 | -5.61058 |
| p=3 | -5.61388 | -5.61292 | -5.61310 | -5.61115 | -5.60915 |
| p=4 | -5.61223 | -5.61091 | -5.61097 | -5.60960 | -5.60764 |
| p=5 | -5.61040 | -5.60892 | -5.60897 | -5.60764 | -5.60605 |
According to the table above, GARCH(1,1) model under t-distributed white noise assumption has the smallest AIC value.
We further use QQ-plot and ACF plot to check whether the residuals of this model satisfy our assumptions.
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")
In the QQ plot, a majority of points lie on the line, with a few of them forming long tails on left and right side. Although the fit in the tails is not absolutely ideal, it is already a significant improvement over the two previous models. In the ACF plot, auto-correlation coefficients at lag 13, 24 and 30 fall on the dividing line of the confidence band, and others are insignificant.
likelihood(nasdaq_garch11_t)## [1] 3538.768log(likelihood(nasdaq_garch11_t))## [1] 8.171534We may conclude that GARCH (1,1) model with t-distributed white noise is a better fit. This model has a likelihood of 3538.768 and a log likelihood of 8.1715. Compared with the result of ARMA(4,4) model and GARCH(4,1) model under normal white noise assumption, this model has a higher likelihood value.
To obtain a better maximum likelihood value, we go a step further to combine ARMA model and GARCH model together. In detail, the only difference between this model and the traditional ARMA model is that \(\{\epsilon_n\}\) is no longer white noise, but a GARCH model.[3]
Using AIC criterion, we select \(p\) and \(q\) values in ARMA(p,q) + GARCH(1,1) model under t distribution assumption.
| q=0 | q=1 | q=2 | q=3 | q=4 | |
|---|---|---|---|---|---|
| p=0 | -5.61679 | -5.61707 | -5.61616 | -5.61660 | -5.61507 |
| p=1 | -5.61695 | -5.61777 | -5.61609 | -5.61528 | -5.61441 |
| p=2 | -5.61585 | -5.61610 | -5.61535 | -5.61653 | -5.61309 |
| p=3 | -5.61645 | -5.61547 | -5.61647 | -5.61527 | -5.62579 |
| p=4 | -5.61535 | -5.61444 | -5.61326 | -5.62611 | -5.62624 |
As illustrated in the table, ARMA(4,4) + GARCH(1,1) has the lowest AIC value.
likelihood(arma44garch11)## [1] 3550.09log(likelihood(arma44garch11))## [1] 8.174728The corresponding likelihood and log likelihood of this model is 3550.09 and 8.1747, which are higher than GARCH(1,1) model.
The leverage \(R_n\) as defined, is the correlation between asset return on day \(n-1\) and the increase in the log volatility from day \(n-1\) to day \(n\). The \(R_n\) can be modeled as: \[ R_n = \frac{e^{2G_n}-1}{e^{2G_n}+1} \] Where \(G_n\) is the usual, Gaussian random walk. Denote the demeaned log return as \(Y_n\), modeling it with: \[ \begin{aligned} Y_n &= e^{\frac{H_n}{2}}\epsilon_n\\ H_n &= \mu_{h}(1-\phi)+\phi H_{n-1}+\beta R_ne^{-\frac{H_{n-1}}{2}}+\omega_{n}\\ G_n &= G_{n-1}+\nu_n \end{aligned} \]
Where \(\beta_n = Y_n\sigma_\eta\sqrt{1-\phi^2}\). \(\epsilon_n\) is an iid \(N(0,1)\) sequence, \(\nu_n\) is an iid \(N(0, \sigma_{\nu}^2)\), and \(\omega_n\) is \(N(0, \sigma_{\omega, n}^2)\) with: \(\sigma^2_{\omega, n}=\sigma^2_{\eta}(1-\phi^2)(1-R_n^2)\)
##
## Attaching package: 'pomp'## The following object is masked from 'package:purrr':
##
## mapFirstly, we need to simulate our model based on an initial parameter setup. As we can see below, the model with such parameter doesn’t fit well to the observed demeaned log return.
params_test <- c(
sigma_nu = exp(-3.5),
mu_h = -0.025,
phi = expit(4),
sigma_eta = exp(-0.07),
G_0 = 0,
H_0=0
)
We will use 3 level parameters to fit our model.
run_level <- 3
NADQ_Np <- switch(run_level, 50, 1e3, 2e3)
NADQ_Nmif <- switch(run_level, 5, 100, 200)
NADQ_Nreps_eval <- switch(run_level, 4, 10, 20)
NADQ_Nreps_local <- switch(run_level, 5, 20, 20)
NADQ_Nreps_global <- switch(run_level, 5, 20, 100)## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: iterators## Loading required package: rngtools## est se
## -1.824468e+03 9.647541e-02NADQ_rw.sd_rp <- 0.02
NADQ_rw.sd_ivp <- 0.1
NADQ_cooling.fraction.50 <- 0.5
NADQ_rw.sd <- rw_sd(sigma_nu = 0.02,
mu_h = 0.02,
phi = 0.02,
sigma_eta = 0.02,
G_0 = ivp(0.1),
H_0 = ivp(0.1)
)summary(r.if1$logLik,digits=5)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3447 3460 3462 3466 3466 3509The maximum likelihood is 3510. Through such local search, the optimal value of \(\sigma_\nu\)is roughly between (0, 0.005), the optimal value of \(\mu_h\) is around -10, the optimal value of \(\phi\)is roughly between (0.8, 0.85) and the optimal value of \(\sigma_eta\) is roughly between (0, 50).
By trying different starting point we will be able to find the global optimized parameters for our model.
pairs(~logLik+sigma_nu+mu_h+phi+sigma_eta,
data=subset(r.if1,logLik>max(logLik)-300))
Create a box to search the best parameters globally.
NADQ_box <- rbind(
sigma_nu=c(0.005,0.05),
mu_h =c(-1,0),
phi = c(0.95,0.99),
sigma_eta = c(0.5,1),
G_0 = c(-2,2),
H_0 = c(-1,1)
)stew(file=paste0("box_eval_",run_level,".rda"),{
if.box <- foreach(i=1:NADQ_Nreps_global,
.packages='pomp',.combine=c, .export=c('if1', 'NADQ_box')) %dopar% {mif2(if1[[1]],
params=apply(NADQ_box,1,function(x)runif(1,x)))}
L.box <- foreach(i=1:NADQ_Nreps_global,
.packages='pomp',.combine=rbind,.export = c('NADQ_Nreps_eval', 'NADQ.filt', 'if', 'NADQ_Np')) %dopar% {
logmeanexp(replicate(NADQ_Nreps_eval, logLik(pfilter(
NADQ.filt,params=coef(if.box[[i]]),Np=NADQ_Np))),
se=TRUE)}
})
timing.box <- .system.time["elapsed"]
r.box <- data.frame(logLik=L.box[,1],logLik_se=L.box[,2],
t(sapply(if.box,coef)))
if(run_level>1) write.table(r.box,file="NADQ_params.csv",
append=TRUE,col.names=FALSE,row.names=FALSE)
summary(r.box$logLik,digits=5)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3446 3459 3463 3469 3471 3513
The iteration shows that the global maximum log likelihood is 3510.
From the figure shown above, we can see that \(\sigma_{\eta}\) converges at near 0.005; \(\sigma_{\nu}\) converges at near 0; \(G_0\) converges at near 0.5; \(\phi\) converges at near 0.8. Meanwhile, \(H_0\) and \(\mu_h\) do not converge.
Some typical time series models, including ARMA, GARCH, and ARMA+GARCH, are first fit to obtain a benchmark likelihood. Using AIC as a criterion for parameter selection, we find that ARMA(4,4) is the best ARIMA model, with a likelihood of 3324.91. Additionally, GARCH(4,1) with normally distributed white noise has the lowest AIC value and has a likelihood of 3476.553. However, both of these two models have heavy-tailed residuals, indicating inappropriate tail fitting. Therefore, GARCH model with t distributed white noise is tried out subsequently, and among various choices of \(p\) and \(q\) values, GARCH(1,1) has the lowest AIC value, generating a likelihood of 3538.77. To further improve the likelihood value, we combine ARMA model and GARCH model together, using ARMA to fit the trend and using GARCH to fit the noise. This yields a larger likelihood of 3550.09, which is obtained by ARMA(4,4) + GARCH(1,1) model.
After this, we use POMP model to fit the NASDAQ data. Through the global searching process, our estimated parameters provide a likelihood of 3510, which performs worse than ARMA+GARCH model, and some of the estimated parameters do not converge. There is definitely something needs to be done to enhance the performance.
We did not consider the seasonality when selecting the appropriate model.
Broader parameters selection needs to be considered when conducting global searching of fitting POMP model.
[1] https://en.wikipedia.org/wiki/Nasdaq_Composite
[2] Analysis of Time Series. Chapter 16: A case study of financial volatility and a POMP model with observations driving latent dynamics. https://ionides.github.io/531w24/16/slides.pdf
[3] David, Ruppert, and S. Matteson David. “Statistics and data analysis for financial engineering: with R examples.” (2015). Section 14.7
[4] David, Ruppert, and S. Matteson David. “Statistics and data analysis for financial engineering: with R examples.” (2015). Section 14.8
[5] Analysis of Time Series. Chapter 16: A case study of financial volatility and a POMP model with observations driving latent dynamics.