---
title: 'stochvolTMB: An R-package for likelihood estimation of stochastic volatility models'
tags:
- R
- maximum likelihood
- laplace approximation
- volatility
- finance
date: "8 December 2020"
output: rticles::joss_article
authors:
- name: Jens Christian Wahl
  orcid: 0000-0002-3812-5111
  affiliation: 1
- name: Norwegian Computing Center
  index: 1
year: 2020
formatted_doi: XX.XXXXX/joss.XXXXX
bibliography: paper.bib
citation_author: Wahl
journal: JOSS
---


# Summary

Stochastic volatility (SV) models are often used to model financial returns that exhibit time-varying and autocorrelated 
variance. The first SV model was introduced by @Taylor1982 and models the logarithm of the variance as a latent autoregressive process of 
order one. Parameter estimation of stochastic volatility models can be challenging and a variety of methods have been 
proposed, such as simulated likelihood [@Liesenfeld2006], quasi-maximum likelihood [@Harvey1994] and 
Markov Chain Monte Carlo methods (MCMC) [@Shepard1998; @Kastner2016]. `stochvolTMB` takes a frequentist approach and estimates
the parameters using maximum likelihood, similar to @Skaug2014. The latent variables are integrated out using the Laplace approximation. 
The models are implemented in `C++` using the `R`-package [@rCore] `TMB` [@TMB2016] for fast and efficient estimation. `TMB` utilizes 
the `Eigen` library [@Eigen2010] for numerical linear algebra and `CppAD` [@CppAD2005] for automatic differentiation of 
the negative log-likelihood. This can lead to substantial speed-up compared to MCMC methods. 

# Implementation

`stochvolTMB` implements stochastic volatility models of the form

\begin{equation}
    \begin{aligned}
        y_t &= \sigma_y e^{h_t/2} \epsilon_t, \quad t = 1, \dots, T, \\
        h_{t+1} &= \phi h_{t} + \sigma_h \eta_t, \quad t = 1, \dots, T-1, \\
        \eta_t &\stackrel{\text{iid}}{\sim} \mathcal{N}(0,1), \\
        \epsilon_t &\stackrel{\text{iid}} {\sim}  F, \\
        h_1 &\sim \mathcal{N} \bigg (0, \frac{\sigma_h}{\sqrt{(1 - \phi^2)}} \bigg )
    \end{aligned}
\end{equation}
where $y_t$ is the observed log return for day $t$, $h_t$ is the logarithm of the conditional variance of day $t$ and $\boldsymbol{\theta} = (\phi, \sigma_y, \sigma_h)$ are the fixed parameters. 
Four distributions are implemented for $\epsilon_t$: (1) The standard normal distribution; (2) The t-distribution with $\nu$ degrees of freedom; 
(3) The skew-normal distribution with skewness parameter $\alpha$; and (4) The leverage model where $(\epsilon_t, \eta_t)$ are both standard normal with correlation parameter
coefficient $\rho$. The last three distributions add an additional fixed parameter to $\boldsymbol{\theta}$. `stochvolTMB` also supports generic functions such as `plot`, `summary`, `predict` and `AIC`. The plotting is 
implemented using `ggplot2` (@ggplot2) and data processing utilizes the `R`-package `data.table` [@datatableRpackage]. 

The parameter estimation is done in an iterative two-step procedure: (1) Optimize the joint negative log-likelihood 
with respect to the latent log-volatility $\boldsymbol{h} = (h_1, \ldots, h_T)$ holding $\boldsymbol{\theta}$ fixed, and (2) Optimizing 
the Laplace approximation of the joint negative log-likelihood w.r.t $\boldsymbol{\theta}$. This procedure is iterated until convergence. 
Standard deviations for the log-volatility and the fixed parameters are obtained by the delta-method [@TMB2016].


<!-- As an example we compare the different models on log-returns for the S&P index from 2005 to 2018: -->

<!-- ```{r warning=FALSE, message=FALSE, } -->
<!-- library(stochvolTMB) -->
<!-- data(spy) -->
<!-- gaussian = estimate_parameters(spy$log_return, model = "gaussian", silent = TRUE) -->
<!-- t_dist = estimate_parameters(spy$log_return, model = "t", silent = TRUE) -->
<!-- skew_gaussian = estimate_parameters(spy$log_return, model = "skew_gaussian", silent = TRUE) -->
<!-- leverage = estimate_parameters(spy$log_return, model = "leverage", silent = TRUE) -->
<!-- ``` -->

<!-- To compare competing models we can use model selection tools such as AIC (@akaike1998): -->

<!-- ```{r message=FALSE} -->
<!-- AIC(gaussian,  -->
<!--     t_dist,  -->
<!--     skew_gaussian,  -->
<!--     leverage) -->
<!-- ``` -->

<!-- Clearly the leverage model is preferred in this example. Notice that the Gaussian model performs the worst and shows the  -->
<!-- importance of having more flexible distributions, even after controlling for the volatility.  -->


`stochvolTMB` differs from `R`-package `stochvol` [@Kastner2016] as `stochvol` performs Bayesian inference using MCMC.
By using optimization instead of simulation we are able to obtain a 5-10 times speed up, depending on the data, model 
and number of observations. 


# References