Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC


Produced with R version 3.4.3 and pomp version 1.16.


14.1 Objectives

  1. Introduce financial volatility and some models used to study it.

  2. Define ARCH, GARCH, and stochastic volatility models.

  3. Provide a case study in using pomp to study a stochastic volatility model.

  4. Discuss this case study as an example of a broader class of nonlinear POMP models, derived from adding stochastically time varying parameters to a linear or nonlinear time series model.

  5. Discuss this case study as an example of an extension of the POMP framework which allows feedback from the observation process to the state process.


14.2 Introduction




14.3 Notation and data

dat <- read.table("sp500.csv",sep=",",header=TRUE)
N <- nrow(dat)
sp500 <- dat$Close[N:1] # data are in reverse order in sp500.csv
par(mfrow=c(1,2))
plot(sp500,type="l")
plot(log(sp500),type="l")




14.3.1 Question. Fitting a stationarity model to stock market returns

  • Look at the plots in Section 3.5.

  • Is it appropriate to fit a stationary model to \({y_{2:N}^*}\), or do we have evidence for violation of stationarity? Explain.




14.4 ARCH and GARCH models

require(tseries)
fit.garch <- garch(sp500,grad = "numerical", trace = FALSE)
L.garch <- logLik(fit.garch)




14.5 Financial leverage



14.5.1 Comment. Time-varying parameters

  • A special case of this model, with the Gaussian random walk having standard deviation zero, is a fixed leverage model.

  • The POMP framework provides a general approach to time-varying parameters. Considering a parameter as a latent, unobserved random process that can progressively change its value over time (following a random walk, or some other stochastic process) leads to a POMP model. The resulting POMP model is usually nonlinear.

  • Many real-world systems are non-stationary and could be investigated using models with time-varying parameters.

  • Some of the midterm projects discovered examples of this.

  • This is one possible thing to explore in your final project.




  • We fit the time-varying leverage model to demeaned daily returns for the S&P 500 index, using the data that were downloaded, detrended and analyzed by Bretó (2014).

  • Here, we just analyze 2002-2012, though it is also interesting to fit to longer intervals, or to compare fits to different intervals.

load(file="sp500-2002-2012.rda")
plot(sp500.ret.demeaned,xlab="business day (2002-2012)",ylab="demeaned S&P 500 return", type="l")

  • We see high volatility around the 2008 financial crisis, and some other episodes.

  • Following the notation and model representation in equation (4) of Bretó (2014), we propose a model, \[ \begin{align} 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, \end{align} \] 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, and \(\{\omega_n\}\) is an iid \(N(0,\sigma_\omega^2)\) sequence.

  • Here, \(H_n\) is the log volatility.


14.6 Building a POMP model

require(pomp)
sp500_statenames <- c("H","G","Y_state")
sp500_rp_names <- c("sigma_nu","mu_h","phi","sigma_eta")
sp500_ivp_names <- c("G_0","H_0")
sp500_paramnames <- c(sp500_rp_names,sp500_ivp_names)
sp500_covarnames <- "covaryt"
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;
 "
sp500_rproc.sim <- paste(rproc1,rproc2.sim)
sp500_rproc.filt <- paste(rproc1,rproc2.filt)
sp500_initializer <- "
  G = G_0;
  H = H_0;
  Y_state = rnorm( 0,exp(H/2) );
"
sp500_rmeasure <- "
   y=Y_state;
"

sp500_dmeasure <- "
   lik=dnorm(y,0,exp(H/2),give_log);
"
sp500_toEstimationScale <- "
  Tsigma_eta = log(sigma_eta);
  Tsigma_nu = log(sigma_nu);
  Tphi = logit(phi);
"

sp500_fromEstimationScale <- "
  Tsigma_eta = exp(sigma_eta);
  Tsigma_nu = exp(sigma_nu);
  Tphi = expit(phi);
"
sp500.filt <- pomp(data=data.frame(y=sp500.ret.demeaned,
                     time=1:length(sp500.ret.demeaned)),
              statenames=sp500_statenames,
              paramnames=sp500_paramnames,
              covarnames=sp500_covarnames,
              times="time",
              t0=0,
              covar=data.frame(covaryt=c(0,sp500.ret.demeaned),
                     time=0:length(sp500.ret.demeaned)),
              tcovar="time",
              rmeasure=Csnippet(sp500_rmeasure),
              dmeasure=Csnippet(sp500_dmeasure),
              rprocess=discrete.time.sim(step.fun=Csnippet(sp500_rproc.filt),delta.t=1),
              initializer=Csnippet(sp500_initializer),
              toEstimationScale=Csnippet(sp500_toEstimationScale), 
              fromEstimationScale=Csnippet(sp500_fromEstimationScale)
)
expit<-function(real){1/(1+exp(-real))}
logit<-function(p.arg){log(p.arg/(1-p.arg))}
params_test <- c(
     sigma_nu = exp(-4.5),  
     mu_h = -0.25,       
     phi = expit(4),     
     sigma_eta = exp(-0.07),
     G_0 = 0,
     H_0=0
  )

sim1.sim <- pomp(sp500.filt, 
               statenames=sp500_statenames,
               paramnames=sp500_paramnames,
               covarnames=sp500_covarnames,
               rprocess=discrete.time.sim(step.fun=Csnippet(sp500_rproc.sim),delta.t=1)
)

sim1.sim <- simulate(sim1.sim,seed=1,params=params_test)
sim1.filt <- pomp(sim1.sim, 
  covar=data.frame(
    covaryt=c(obs(sim1.sim),NA),
    time=c(timezero(sim1.sim),time(sim1.sim))),
  tcovar="time",
  statenames=sp500_statenames,
  paramnames=sp500_paramnames,
  covarnames=sp500_covarnames,
  rprocess=discrete.time.sim(step.fun=Csnippet(sp500_rproc.filt),delta.t=1)
)




14.7 Filtering on simulated data

run_level <- 3 
sp500_Np <-          c(100,1e3,2e3)
sp500_Nmif <-        c(10, 100,200)
sp500_Nreps_eval <-  c(4,  10,  20)
sp500_Nreps_local <- c(10, 20, 20)
sp500_Nreps_global <-c(10, 20, 100)
require(doParallel)
registerDoParallel()
stew(file=sprintf("pf1.rda",run_level),{
  t.pf1 <- system.time(
    pf1 <- foreach(i=1:sp500_Nreps_eval[run_level],.packages='pomp',
                   .options.multicore=list(set.seed=TRUE)) %dopar% try(
                     pfilter(sim1.filt,Np=sp500_Np[run_level])
                   )
  )
},seed=493536993,kind="L'Ecuyer")
(L.pf1 <- logmeanexp(sapply(pf1,logLik),se=TRUE))
##                          se 
## -3.658756e+03  9.708603e-02




14.8 Fitting the stochastic leverage model to S&P500 data

sp500_rw.sd_rp <- 0.02
sp500_rw.sd_ivp <- 0.1
sp500_cooling.fraction.50 <- 0.5

stew("mif1.rda",{
   t.if1 <- system.time({
   if1 <- foreach(i=1:sp500_Nreps_local[run_level],
                  .packages='pomp', .combine=c,
                  .options.multicore=list(set.seed=TRUE)) %dopar% try(
                    mif2(sp500.filt,
                         start=params_test,
                         Np=sp500_Np[run_level],
                         Nmif=sp500_Nmif[run_level],
                         cooling.type="geometric",
                         cooling.fraction.50=sp500_cooling.fraction.50,
                         transform=TRUE,
                         rw.sd = rw.sd(
                            sigma_nu  = sp500_rw.sd_rp,
                            mu_h      = sp500_rw.sd_rp,
                            phi       = sp500_rw.sd_rp,
                            sigma_eta = sp500_rw.sd_rp,
                            G_0       = ivp(sp500_rw.sd_ivp),
                            H_0       = ivp(sp500_rw.sd_ivp)
                         )
                    )
                  )
    
    L.if1 <- foreach(i=1:sp500_Nreps_local[run_level],.packages='pomp',
                      .combine=rbind,.options.multicore=list(set.seed=TRUE)) %dopar% 
                      {
                        logmeanexp(
                          replicate(sp500_Nreps_eval[run_level],
                                    logLik(pfilter(sp500.filt,params=coef(if1[[i]]),Np=sp500_Np[run_level]))
                          ),
                          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="sp500_params.csv",append=TRUE,col.names=FALSE,row.names=FALSE)
summary(r.if1$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -3990   -3957   -3956   -3959   -3955   -3954
pairs(~logLik+sigma_nu+mu_h+phi+sigma_eta,data=subset(r.if1,logLik>max(logLik)-20))




14.9 Likelihood maximization using randomized starting values.

sp500_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="box_eval.rda",{
  t.box <- system.time({
    if.box <- foreach(i=1:sp500_Nreps_global[run_level],.packages='pomp',.combine=c,
                  .options.multicore=list(set.seed=TRUE)) %dopar%  
      mif2(
        if1[[1]],
        start=apply(sp500_box,1,function(x)runif(1,x))
      )
    
    L.box <- foreach(i=1:sp500_Nreps_global[run_level],.packages='pomp',.combine=rbind,
                      .options.multicore=list(set.seed=TRUE)) %dopar% {
                        set.seed(87932+i)
                        logmeanexp(
                          replicate(sp500_Nreps_eval[run_level],
                                    logLik(pfilter(sp500.filt,params=coef(if.box[[i]]),Np=sp500_Np[run_level]))
                          ), 
                          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="sp500_params.csv",append=TRUE,col.names=FALSE,row.names=FALSE)
summary(r.box$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -3996   -3986   -3957   -3966   -3956   -3954
pairs(~logLik+log(sigma_nu)+mu_h+phi+sigma_eta+H_0,data=subset(r.box,logLik>max(logLik)-10))




14.10 Benchmark likelihoods for alternative models

require(tseries)
fit.garch.benchmark <- garch(sp500.ret.demeaned,grad = "numerical", trace = FALSE)
L.garch.benchmark <- logLik(fit.garch.benchmark)




14.11 Appendix: Proper weighting for a partially plug-and-play algorithm with a perfectly observed state space component.





Acknowledgment

These notes draw on material in a tutorial on financial leverage analysis in pomp.


References

Bretó, C. 2014. On idiosyncratic stochasticity of financial leverage effects. Statistics & Probability Letters 91:20–26.

Cowpertwait, P. S., and A. V. Metcalfe. 2009. Introductory time series with R. Springer Science & Business Media.

Ionides, E. L., D. Nguyen, Y. Atchadé, S. Stoev, and A. A. King. 2015. Inference for dynamic and latent variable models via iterated, perturbed Bayes maps. Proceedings of the National Academy of Sciences of the U.S.A. 112:719–724.

Liu, J. S., and R. Chen. 1998. Sequential Monte Carlo methods for dynamic systems. Journal of the American Statistical Association 93:1032–1044.