Stochastic volatility of Apple stock

1. Introduction
2. Data Preparing
3. GARCH Modeling
4. POMP model of Carles Breto
5. Conclusions
6. Reference

1. Introduction

Stochastic volatility models in statistics are those in which the variance of a stochastic process is itself randomly distributed\(^{[1]}\). As we all know, returns on investments in large companies are often found to be approximately uncorrelated. The variability of the returns can fluctuate considerably. Understanding this volatility is important for quantifying and managing the risk of investments.\(^{[2]}\)

In this project, I will use stochastic volatility models to analyze Apple stock. Apple is famous for its size and revenues, therefore, to some extent, analyzing Apple stock can represent the stochastic volatility of stocks of large companies. I used two models, the first one is the Garch model, and the other is the POMP model proposed by Carles Breto (2014). Then I compared the two models based on the likelihood values.

2. Data Preparing

In the project, we focused on the Apple stock adjusted close price from 1983 to 2018. The data we used is monthly data. We got the data from https://finance.yahoo.com/. The plots below show the raw data and log-transformed data. From the plots we can see that the Apple stock price increases rapidly, especially for the last ten years.

rm(list = ls())
setwd(file.path(getwd()) )
par(mfrow=c(1,2)) 
data <- read.table("AAPL.csv",sep=",",header=TRUE)
data$Date=as.Date(as.character(data$Date),"%m/%d/%Y")
plot(data$Date, data$AdjClose, xlab="date",ylab="Adj.Close",type="l", main='Apple adjusted close price')
plot(data$Date, log(data$AdjClose),xlab="date",ylab="log Adj.Close",type="l", main='Apple log adjusted close price')

The return refers to the difference of the log price: \[y_n = log(z_n) - log(z_{n-1}) \]

In order to get a mean stationary dataset, we removed the mean of the returns and we will fit models to the demeaned returns. The demeaned returns are shown as follows. We noticed that the volatility become large significantly around 2001.

return = diff(log(data$AdjClose))
demeanReturn = return - mean(return)
plot(data$Date[-1], demeanReturn,type='l', lwd=2, main='Demeaned returns', xlab = 'date', ylab = 'demeaned returns')
abline(h=mean(demeanReturn),col='red', lwd=2)

If we plot the Auto-correlation function plots, we found that the demeaned returns are uncorrelated and the squared demeaned returns have a relatively stronger correlation even though the correlation is still not very strong.

par(mfrow=c(1,2)) 
acf(demeanReturn, main='Demeaned Return')
acf(demeanReturn^2,main='Squared-demeaned Return')

3. GARCH Modeling

The GARCH(p,q) model is often preferred by financial modeling. It has the form:

\[Y_n = \epsilon_n \sqrt{V_n} \]

\[V_n = \alpha_0 + \sum_{j=1}^p \alpha_j Y_{n-j}^2 + \sum_{k=1}^q \beta_k V_{n-k} \]

Where \(\epsilon_{1:N}\) is white noise.

In practice, a popular choice is the GARCH(1.1) model\(^{[3]}\). I used the garchFit function in R to fit the demeaned returns of the Apple stock. In the following plot, the orange line shows the confidence interval of the returns. Here, the confidence interval means \(0\pm qnorm(0.975) * \sqrt{volatility}\). The plot shows that most of the raw data are listed between the confidence interval generated by Garch(1,1) model.

The log likelihood value of the Garch(1,1) model is 289.6891, we need to compare it with likelihood values of POMP model.

library(knitr)
#require(tseries)
library(fGarch)

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

garch_model <- garchFit(~garch(1,1),demeanReturn,trace=F)

summary(garch_model)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = demeanReturn, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7fa2e2c23cf8>
##  [data = demeanReturn]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1  
## -2.7867e-17   6.4466e-04   1.4124e-01   8.3137e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -2.787e-17   5.416e-03    0.000 1.000000    
## omega   6.447e-04   3.944e-04    1.635 0.102140    
## alpha1  1.412e-01   4.029e-02    3.505 0.000456 ***
## beta1   8.314e-01   4.635e-02   17.938  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  289.6891    normalized:  0.6721326 
## 
## Description:
##  Sat May  2 19:18:34 2020 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  100.6042  0           
##  Shapiro-Wilk Test  R    W      0.9722085 2.607449e-07
##  Ljung-Box Test     R    Q(10)  7.503535  0.6772051   
##  Ljung-Box Test     R    Q(15)  13.03842  0.5993299   
##  Ljung-Box Test     R    Q(20)  16.22327  0.7026802   
##  Ljung-Box Test     R^2  Q(10)  8.058676  0.6231055   
##  Ljung-Box Test     R^2  Q(15)  9.233138  0.8650152   
##  Ljung-Box Test     R^2  Q(20)  11.87188  0.9204152   
##  LM Arch Test       R    TR^2   9.371011  0.6709568   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -1.325704 -1.287967 -1.325874 -1.310804

sd=garch_model@sigma.t

plot(as.Date(data$Date)[-1], demeanReturn, ylab="Returns", xlab = 'Date', type = 'l', main = 'Garch(1,1) and fitted CI')
lines(as.Date(data$Date)[-1], qnorm(0.025)*sd, lty=2, col='orange')
lines(as.Date(data$Date)[-1], qnorm(0.975)*sd, lty=2, col='orange')
legend('bottomright', c('returns','Confidence intervals'), col = c('black','orange'), lty = c(1,2))

Next, I did some residual checking. The ACF plot shows the residuals are almost uncorrelated but have relative strong correlations at lag 10, 13, and 21. The QQ plot shows that the residuals are nearly normal distributed but left-skewed. These are indicators of the not perfect fitting.

par(mfrow=c(1,2)) 
acf(garch_model@residuals)
qqnorm(garch_model@residuals)
qqline(garch_model@residuals,col=2)

4. POMP model of Carles Breto

4.1 Financial leverage

The leverage \(R_n\) on day n is the correlation between return on day n-1 and the increase in the log volatility from day n-1 to day n.^{[2]}

Define: \[ R_n = \frac{exp\{2G_n\} -1 } {exp\{2G_n\} + 1} \]

Where \(\{ G_n\}\) is the Gaussian random walk.

4.2 Model description

The POMP model allows the parameters varying over time. The parameters are latent random processes.

Carles Breto(2014) proposed a model that can be constructed by POMP framework:

\[Y_n = exp \{ H_n /2\} \epsilon_n \]

\[H_n = \mu_h (1-\phi) + \phi H_{n-1} + \beta_{n-1} R_n exp\{-H_{n-1} /2 \} + \omega_n \] \[G_n = G_{n-1} + \nu_n \] 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)\) sequence, \(\{\omega_n \}\) is an iid \(N(0, \sigma_{\omega}^2)\) sequence.

\[\sigma_{\omega}^2 = \sigma_{\eta}^2 (1-\phi^2)(1-R_n^2) \] Here, \(H_n\) is the log volatility.

There are five parameters: \(\mu_h\), \(\phi\), \(\sigma_{\eta}\), \(\sigma_\nu\), and \(R_n\)

4.3 building POMP model

In this model, the state variable is defined as: \[X_n = (G_n , H_n, Y_n)\]

The measurement process is the perfect observation of the \(Y_n\) component of the state space.

For iterated filtering, it is convenient to transform parameters to be defined on the whole real line, so we transform \(\sigma_\eta\) into \(log(\sigma_\eta)\), transform \(\sigma_\nu\) into \(log(\sigma_\nu)\), and transform \(\phi\) into \(logit(\phi)\) since \(\phi \in [0,1]\)

set.seed(594709947L)
library(ggplot2)
theme_set(theme_bw())
library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ✓ purrr   0.3.4

## Warning: package 'tibble' was built under R version 3.6.2

## Warning: package 'purrr' was built under R version 3.6.2

## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks timeSeries::filter(), stats::filter()
## x dplyr::lag()    masks timeSeries::lag(), stats::lag()

library(plyr)

## -------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## -------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:purrr':
## 
##     compact

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

library(foreach)

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

#library(doMC)
library(pomp)

## Warning: package 'pomp' was built under R version 3.6.2

## Welcome to pomp version 2.7.1.0!
## As of version 2.7.1.0, important changes have been made to the
## default settings of the particle filtering algorithms in
## 'pfilter', 'mif2', 'pmcmc', 'bsmc2'.
## These changes are not backward compatible.
## See the package NEWS for the details.
## 
## For information on upgrading your pomp version < 2 code, see the
## "pomp version 2 upgrade guide" at https://kingaa.github.io/pomp/.

## 
## Attaching package: 'pomp'

## The following object is masked from 'package:purrr':
## 
##     map

## The following object is masked from 'package:timeSeries':
## 
##     time<-

stopifnot(packageVersion("pomp")>="2.0")

#names
apple_statenames <- c("H","G","Y_state")
apple_rp_names <- c("sigma_nu","mu_h","phi","sigma_eta")
apple_ivp_names <- c("G_0","H_0")
apple_paramnames <- c(apple_rp_names,apple_ivp_names)

#rprocess
rproc1 <- "
  double beta,omega,nu;
  omega = rnorm(0,sigma_eta * sqrt( 1- phi*phi ) * 
    sqrt(1-tanh(G)*tanh(G)));
  nu = rnorm(0, sigma_nu);
  G += nu;
  beta = Y_state * sigma_eta * sqrt( 1- phi*phi );
  H = mu_h*(1 - phi) + phi*H + beta * tanh( G ) 
    * exp(-H/2) + omega;
"
rproc2.sim <- "Y_state = rnorm( 0,exp(H/2) );"

rproc2.filt <- "Y_state = covaryt;"
apple_rproc.sim <- paste(rproc1,rproc2.sim)
apple_rproc.filt <- paste(rproc1,rproc2.filt)


#rinit
apple_rinit <- "
  G = G_0;
  H = H_0;
  Y_state = rnorm( 0,exp(H/2) );"

#measure
apple_rmeasure <- "y=Y_state;"
apple_dmeasure <- " lik=dnorm(y,0,exp(H/2),give_log);"

# transforms
apple_partrans <- parameter_trans(
  log=c("sigma_eta","sigma_nu"),
  logit="phi")

apple.filt <- pomp(data=data.frame(
    y=demeanReturn,time=1:length(demeanReturn)),
  statenames=apple_statenames,
  paramnames=apple_paramnames,
  times="time",
  t0=0,
  covar=covariate_table(
    time=0:length(demeanReturn),
    covaryt=c(0,demeanReturn),
    times="time"),
  rmeasure=Csnippet(apple_rmeasure),
  dmeasure=Csnippet(apple_dmeasure),
  rprocess=discrete_time(step.fun=Csnippet(apple_rproc.filt), delta.t=1),
  rinit=Csnippet(apple_rinit),
  partrans=apple_partrans
)

To test and develop the code, we first simulated from the model.

## sim_pomp
params_test <- c(
  sigma_nu = exp(-4.5),  
  mu_h = -5,     
  phi = expit(-3),   
  sigma_eta = exp(-0.07),
  G_0 = 0,
  H_0=0
)
  
sim1.sim <- pomp(apple.filt, 
  statenames=apple_statenames,
  paramnames=apple_paramnames,
  rprocess=discrete_time(
    step.fun=Csnippet(apple_rproc.sim),delta.t=1)
)

sim1.sim <- simulate(sim1.sim,seed=1,params=params_test)

plot(Y_state~time, data=sim1.sim, type='l', col='blue', main="Observed returns and simulated returns", ylab="Returns")
lines(demeanReturn,col='black')
legend('topright' , c("Observed Returns","Simulated Returns"), col=c("black","blue"), lty=c(1,1))

The simulated returns and observed returns are shown above. In general, the simulated data coincides with the volatility of the observed data.

4.4 Local search for maximizing the likelihood for Apple stock data

Now, I use the IF2 algorithm of Ionides et al. (2015) to fit the stochastic leverage model to Apple stock data. We can get the log likelihood value as well as convergence diagnostics plots.

library(doParallel)

## Loading required package: iterators

## Loading required package: parallel

registerDoParallel(6)
library(doRNG)

## Loading required package: rngtools

registerDoRNG(34118892)

run_level <- 3 
apple_Np <-           switch(run_level, 100, 1e3, 2e3)
apple_Nmif <-         switch(run_level,  10, 100, 200)
apple_Nreps_eval <-   switch(run_level,   4,  10,  20)
apple_Nreps_local <-  switch(run_level,  10,  20,  20)
apple_Nreps_global <- switch(run_level,  10,  20, 100)

## mif_setup 
apple_rw.sd_rp <- 0.02
apple_rw.sd_ivp <- 0.1
apple_cooling.fraction.50 <- 0.5
apple_rw.sd <- rw.sd(
  sigma_nu  = apple_rw.sd_rp,
  mu_h      = apple_rw.sd_rp,
  phi       = apple_rw.sd_rp,
  sigma_eta = apple_rw.sd_rp,
  G_0       = ivp(apple_rw.sd_ivp),
  H_0       = ivp(apple_rw.sd_ivp)
)    


## mif
stew(file=sprintf("mif1-%d.rda",run_level),{
  t.if1 <- system.time({
  if1 <- foreach(i=1:apple_Nreps_local,
    .packages='pomp', .combine=c) %dopar% mif2(apple.filt,
      params=params_test,
      Np=apple_Np,
      Nmif=apple_Nmif,
      cooling.fraction.50=apple_cooling.fraction.50,
      rw.sd = apple_rw.sd)
  L.if1 <- foreach(i=1:apple_Nreps_local,
    .packages='pomp', .combine=rbind) %dopar% logmeanexp(
      replicate(apple_Nreps_eval, logLik(pfilter(apple.filt,
        params=coef(if1[[i]]),Np=apple_Np))), se=TRUE)
  })
},seed=318817883,kind="L'Ecuyer")


r.if1 <- data.frame(logLik=L.if1[,1],logLik_se=L.if1[,2],
  t(sapply(if1,coef)))
if (run_level>1) write.table(r.if1,file="apple_params.csv",
  append=TRUE,col.names=FALSE,row.names=FALSE)

summary(r.if1$logLik,digits=5)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   298.1   300.9   301.1   300.7   301.2   301.4

The maximum log likelihood value of this POMP model is 301.4, which is a little bigger than 289.6891, the log likelihood of Garch model. However, the POMP model has more parameters. According to the AIC criteria \(AIC = 2k - 2 log(\hat{L})\) the two models actually perform similarly.

pairs(~logLik+sigma_nu+mu_h+phi+sigma_eta, data=subset(r.if1,logLik>max(logLik)-20))

The plot above shows us the geometry of the likelihood surface in a neighborhood of this point estimate. Because we only use fixed starting values, the points seems sparse.

plot(if1)

Next, I checked the convergence diagnostics of the POMP model. From the MIF2 convergence diagnostics, the log likelihood converges after 50 iterations. Besides, the \(\sigma_{\nu}\) seems to converge after 150 iterations. The parameter \(\phi\) gathers into a small interval. The \(G_0\) and \(H_0\) seems to converge after 100 iterations. But the parameters \(\mu_h\) and \(\sigma_\eta\) seem not converge.

In the next section, I will use randomized starting values to implement the POMP model again.

4.5 Global search for maximizing the likelihood

Using randomly starting values from a large box is useful for global maximization. I chose a box that might contain the plausible parameter values. The box is shown as follows.

## box 
apple_box <- rbind(
 sigma_nu=c(0.0,0.01),
 mu_h    =c(-5,-3),
 phi = c(0.93,0.98),
 sigma_eta = c(0,2),
 G_0 = c(-2,2),
 H_0 = c(-1,1)
)

## box_eval
stew(file=sprintf("box_eval-%d.rda",run_level),{
  t.box <- system.time({
    if.box <- foreach(i=1:apple_Nreps_global,
      .packages='pomp',.combine=c) %dopar% mif2(if1[[1]],
        params=apply(apple_box,1,function(x)runif(1,x)))
    L.box <- foreach(i=1:apple_Nreps_global,
      .packages='pomp',.combine=rbind) %dopar% {
         logmeanexp(replicate(apple_Nreps_eval, logLik(pfilter(
         apple.filt,params=coef(if.box[[i]]),Np=apple_Np))), 
           se=TRUE)
       }
  })
},seed=290860873,kind="L'Ecuyer")

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="apple_params.csv",
  append=TRUE,col.names=FALSE,row.names=FALSE)
summary(r.box$logLik,digits=5)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   293.2   298.2   300.2   299.7   301.1   301.4

I noticed that the best log likelihood value is still 301.4, which is the same with the log likelihood of the POMP model with given starting parameter values. This gave us some confidence that the MLE is reasonable.

pairs(~logLik+log(sigma_nu)+mu_h+logit(phi)+log(sigma_eta), data=subset(r.box,logLik>max(logLik)-10))

There are more estimated points when we use random starting parameter values. This plot gives us the feeling of the global geometry of the likelihood surface. To reach the maximum of the log likelihood, the best \(log(\sigma_\nu)\) fall in the range (-9,-5), the best \(\mu_h\) is around -4, the best \(log(\phi)\) is around 3, the best \(log(\sigma_\eta)\) is around 0.

plot(if.box)

The convergence diagnostics is similar to the convergence diagnostics in the last section. It is obvious that the loglikelihood converges after 50 iterations. The \(\sigma_\mu\) and \(\phi\) approximately converge to 0 and 0.95 respectively. \(\mu_h\) and \(\sigma_\eta\) seem not to converge.

5. Conclusions

In the project, we used two stochastic volatility models, Garch model and POMP model, to fit the Apple stock. I chose the most popular Garch(1,1) model to fit the data and got the log likelihood value 289.7. The POMP model with fixed or random starting parameters values have log likelihood value 301.4. Since the Garch(1,1) model has 3 fixed parameters while the POMP model has 6 fixed parameters, AIC criteria favors the POMP model.

The convergence diagnostics of POMP model show that the log likelihood converges well. This is a good indicator that the POMP model is not bad for the Apple Stock data. There are two parameters in POMP model that seem not to converge: \(\mu_h\) and \(\sigma_\eta\). Maybe we need to fine-tune the starting parameter box. Or maybe we need to iterate more times.