Estimation of realized stochastic volatility models using Hamiltonian Monte Carlo-Based methods

Abstract

This study develops and compares performance of Hamiltonian Monte Carlo (HMC) and Riemann manifold Hamiltonian Monte Carlo (RMHMC) samplers with that of multi-move Metropolis-Hastings sampler to estimate stochastic volatility (SV) and realized SV models with asymmetry effect. In terms of inefficiency factor, empirical results show that the RMHMC sampler give the best performance for estimating parameters, followed by multi-move Metropolis-Hastings sampler. In particular, it is also shown that RMHMC sampler offers a greater advantage in the mixing property of latent volatility chains and in the computational time than HMC sampler.

Appendices

Appendix 1: Estimation of RASV model

Let us take the logarithm of the full conditional posterior density defined in Eq. (3):

$$\begin{aligned} \ln p(\theta _1,\varvec{\alpha }|\mathbf {R},\mathbf {Y})&= constant +\sum _{t=1}^{T}\left( -\frac{1}{2}\alpha _t -\frac{1}{2}\widetilde{R}_t^2 \right) -T\ln (\sigma _\epsilon \sigma _u) \nonumber \\&\quad -\frac{1}{2\sigma _u^2}\sum _{t=1}^{T}\left( Y_t-c-\alpha _t \right) ^2 -T\ln (\sigma _\eta ) + \frac{1}{2}\ln (1-\phi ^2) \nonumber \\&\quad -\frac{T-1}{2}\ln (1-\rho ^2) -\frac{1-\phi ^2}{2\sigma _\eta ^2}\alpha _1^2 -\frac{1}{2\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t^2 \nonumber \\&\quad -\frac{1}{2}\ln (v_c) -\frac{(c-m_c)^2}{2v_c} -\left( \frac{a_{\sigma _u}}{2}+1\right) \ln \left( \sigma _u^2\right) -\frac{b_{\sigma _u}}{2\sigma _u^2} \nonumber \\&\quad -\left( \frac{a_{\sigma _\epsilon }}{2}+1\right) \ln \left( \sigma _\epsilon ^2\right) -\frac{b_{\sigma _\epsilon }}{2\sigma _\epsilon ^2} +(a_\phi -1)\ln (1+\phi ) \nonumber \\&\quad +(b_\phi -1)\ln (1-\phi ) -\left( \frac{a_{\sigma _\eta }}{2}+1\right) \ln \left( \sigma _\eta ^2\right) -\frac{b_{\sigma _\eta }}{2\sigma _\eta ^2} \nonumber \\&\quad +(a_\rho -1)\ln (1+\rho ) +(b_\rho -1)\ln (1-\rho ) \end{aligned}$$

(5)

where

$$\begin{aligned} \widetilde{R}_t&= R_t\sigma _\epsilon ^{-1} e^{-\frac{1}{2}h_t},~~~\hbox {for}~t=1,\ldots ,T,~\hbox {and} \\ \widetilde{\alpha }_t&= \alpha _t-\phi \alpha _{t-1}-\rho \sigma _\eta \widetilde{R}_{t-1},~~~\hbox {for}~t=2,\ldots ,T. \end{aligned}$$

The following subsections describe how to sample latent variable $\varvec{\alpha }$ and parameters in $\theta _1=(c,\sigma _u,\sigma _\epsilon ,\phi ,\sigma _\eta ,\rho )$.

1.1 Updating latent variable $\varvec{\alpha }$

Sampling from the conditional posterior density $p(\varvec{\alpha }|\mathbf {R},\mathbf {Y},\theta _1)$, is a difficult work because of the high-dimensional log-normal structure. Therefore, we sample latent variable $\varvec{\alpha }$ using the HMC and RMHMC algorithms, where the expressions required by both algorithms to sample the entire latent volatility at once are evaluated as follows.

On the basis of Eq. (5), the logarithm of the full conditional posterior density of $\varvec{\alpha }$ has the following expression:

$$\begin{aligned} {\mathcal {L}}(\varvec{\alpha })&= constant +\sum \nolimits _{t=1}^{T}\left( -\frac{1}{2}\alpha _t -\frac{1}{2}\widetilde{R}_t^2 \right) -\frac{1}{2\sigma _u^2}\sum _{t=1}^{T}\left( Y_t-c-\alpha _t \right) ^2\\&\quad -\frac{1-\phi ^2}{2\sigma _h^2}\alpha _1^2 -\frac{1}{2\sigma _h^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t^2. \end{aligned}$$

The important aspect of the HMC and RMHMC algorithms is that both algorithms require the calculation of a gradient vector $\nabla _{\alpha _t}{\mathcal {L}}(\alpha _t)$. These are expressed by

$$\begin{aligned} \nabla _{\alpha _1}{\mathcal {L}}(\varvec{\alpha })&= \nabla _{\alpha _1}{\mathcal {L}}_{\mathbf {R},\mathbf {Y}} -\frac{1-\phi ^2}{\sigma _\eta ^2}\alpha _1 +\frac{1}{\sigma _\eta ^2(1-\rho ^2)}\widetilde{\alpha }_2\left( \phi -\frac{1}{2}\rho \sigma _\eta \widetilde{R}_1\right) , \\ \nabla _{\alpha _T}{\mathcal {L}}(\varvec{\alpha })&= \nabla _{\alpha _T}{\mathcal {L}}_{\mathbf {R},\mathbf {Y}} -\frac{1}{\sigma _\eta ^2(1-\rho ^2)}\widetilde{\alpha }_T, \\ \nabla _{\alpha _i}{\mathcal {L}}(\varvec{\alpha })&= \nabla _{\alpha _i}{\mathcal {L}}_{\mathbf {R},\mathbf {Y}} -\frac{1}{\sigma _\eta ^2(1-\rho ^2)}\left[ \widetilde{\alpha }_i -\widetilde{\alpha }_{i+1}\left( \phi -\frac{1}{2}\rho \sigma _\eta \widetilde{R}_i\right) \right] , \end{aligned}$$

for $i=2,\ldots ,T-1$, and

$$\begin{aligned} \nabla _{\alpha _t}{\mathcal {L}}_{\mathbf {R},\mathbf {Y}} = -\frac{1}{2}+\frac{1}{2} \widetilde{R}_t^2 +\frac{1}{\sigma _u^2}(Y_t-c-\alpha _t), ~~~t=1,\ldots ,T \end{aligned}$$

To apply the leapfrog-RMHMC algorithm, then we also need expressions for the metric tensor and its partial derivatives. The metric tensor $\mathbf {M}(\varvec{\alpha })$ is a symmetric tridiagonal matrix, where elements of main diagonal are calculated as follows:

$$\begin{aligned} \mathbf {M}(1,1)&= \frac{1}{2} +\frac{1}{\sigma _u^2} +\frac{1-\phi ^2}{\sigma _\eta ^2} +\frac{1}{\sigma _\eta ^2(1-\rho ^2)}\left( \phi ^2+\frac{1}{4}\rho ^2\sigma _\eta ^2\right) ,\\ \mathbf {M}(i,i)&= \frac{1}{2} +\frac{1}{\sigma _u^2} +\frac{1}{\sigma _\eta ^2(1-\rho ^2)}\left( 1+\phi ^2+\frac{1}{4}\rho ^2\sigma _\eta ^2\right) ,\\ \mathbf {M}(T,T)&= \frac{1}{2} +\frac{1}{\sigma _u^2} +\frac{1}{\sigma _\eta ^2(1-\rho ^2)} \end{aligned}$$

for $i=2,\ldots ,T-1$, and elements of superdiagonal and subdiagonal are calculated as

$$\begin{aligned} \mathbf {M}(i-1,i) = \mathbf {M}(i,i-1) = -\frac{\phi }{\sigma _\eta ^2(1-\rho ^2)}, \end{aligned}$$

for $i=2,\ldots ,T$. Since the above metric tensor is fixed, that is not a function of the latent volatility $\varvec{\alpha }$, the associated partial derivatives with respect to the latent volatility are zero. This implies that both metric tensor and its inverse are not recalculated in the generalized leapfrog integration scheme.

1.2 Updating parameter $\phi$

As far as $\phi$ is only concerned to be sampled, the logarithm of conditional posterior density for $\phi$ on the basis of Eq. (5) has the following expression:

$$\begin{aligned} {\mathcal {L}}(\phi )&= constant + \frac{1}{2}\ln (1-\phi ^2) -\frac{\alpha _1^2}{2\sigma _\eta ^2}(1-\phi ^{2}) -\frac{1}{2\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t^2\\&\quad +\,(a_\phi -1)\ln (1+\phi ) +(b_\phi -1)\ln (1-\phi ), \end{aligned}$$

which is in non-standard form, and therefore we cannot sample from it directly. In the following, we evaluate the expressions required by the HMC and RMHMC algorithms to sample $\phi$.

For dealing with the constraints $-1<\phi <1$, it is necessary to implement the transformations $\phi =\tanh (\tilde{\phi })$ and then introduce a Jacobian factor into the acceptance ratio. The derivative of the transformation is $\hbox {d}\phi /\hbox {d}\tilde{\phi }=1-\tanh ^2(\tilde{\phi })=1-\phi ^2$. Then, the partial derivative of the above log joint posterior with respect to the transformed parameter is as follows:

$$\begin{aligned} \nabla _{\tilde{\phi }}{\mathcal {L}}(\phi )&= -\phi + \frac{\alpha _1^2}{\sigma _\eta ^2}(\phi -\phi ^3) +\frac{1-\phi ^2}{\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T} \widetilde{\alpha }_t\alpha _{t-1} \\&\quad +\,(a_\phi -1)(1-\phi ) -(b_\phi -1)(1+\phi ). \end{aligned}$$

To apply the leapfrog RMHMC algorithm, then we also need expressions for the metric tensor and its partial derivative with respect to the transformed parameter. The metric tensor $\mathbf {M}(\phi )$ for the log joint posterior density is calculated as follows:

$$\begin{aligned} \mathbf {M}(\phi ) = 2\phi ^2 +\frac{(T-1)(1-\phi ^2)}{1-\rho ^2} +(a_\phi +b_\phi -2)(1-\phi ^2). \end{aligned}$$

Thus the partial derivative of metric tensor $\mathbf {M}(\phi )$ follows as

$$\begin{aligned} \frac{\partial \mathbf {M}(\phi )}{\partial \tilde{\phi }}&= 2\phi (1-\phi ^2)\left( 4-\frac{T-1}{1-\rho ^2}-a_\phi -b_\phi \right) . \end{aligned}$$

1.3 Updating parameters $(\sigma _\epsilon ,\sigma _\eta ,\rho )$

On the basis of Eq. (5), the log joint posterior density of ${\varSigma } = (\sigma _\epsilon ,\sigma _\eta ,\rho )$ has the following expression:

$$\begin{aligned} {\mathcal {L}}({\varSigma })&= constant-T\ln (\sigma _\epsilon \sigma _\eta ) -\frac{1}{2}\sum _{t=1}^{T}\widetilde{R}_t^2-\frac{(1-\phi ^2)\alpha _1^2}{2\sigma _\eta ^2} \\&\quad -\frac{T-1}{2}\ln (1-\rho ^2) -\frac{1}{2\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t^2\\&\quad -\left( a_{\sigma _\epsilon }+2\right) \ln \left( \sigma _\epsilon \right) -\frac{b_{\sigma _\epsilon }}{2\sigma _\epsilon ^2} -\left( a_{\sigma _\eta }+2\right) \ln \left( \sigma _\eta \right) -\frac{b_{\sigma _\eta }}{2\sigma _\eta ^2} \\&\quad +\,(a_\rho -1)\ln (1+\rho ) +(b_\rho -1)\ln (1-\rho ), \end{aligned}$$

which is in non-standard form. The expressions required by the HMC and RMHMC algorithms to sample ${\varSigma }$ are evaluated as follows.

For dealing with the constraints $-1<\rho <1$ and $\sigma _\epsilon ,\sigma _\eta >0$, it is necessary to implement the transformations $\rho =\tanh (\tilde{\rho })$, $\sigma _\epsilon =\exp (\tilde{\sigma }_\epsilon )$, and $\sigma _\eta =\exp (\tilde{\sigma }_\eta )$ and then introduce a Jacobian factor into the acceptance ratio. Notice that Girolami and Calderhead (2011) did not take a transformation for the parameter $\sigma _\epsilon$. The derivatives of the transformations are $\hbox {d}\rho /\hbox {d}\tilde{\rho }=1-\rho ^2$, $\hbox {d}\sigma _\epsilon /\hbox {d}\tilde{\sigma }_\epsilon =\sigma _\epsilon$, and $\hbox {d}\sigma _\eta /\hbox {d}\tilde{\sigma }_\eta =\sigma _\eta$. Then, the partial derivatives of the above log joint posterior with respect to the transformed parameters are as follows:

$$\begin{aligned} \nabla _{\tilde{\sigma }_\epsilon }{\mathcal {L}}({\varSigma })&= -T +\sum _{t=1}^{T}\widetilde{R}_t^2 -\frac{\rho }{\sigma _\eta (1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t\widetilde{R}_{t-1} -\left( a_{\sigma _\epsilon }+2\right) +\frac{b_{\sigma _\epsilon }}{\sigma _\epsilon ^2} \\ \nabla _{\tilde{\sigma }_\eta }{\mathcal {L}}({\varSigma })&= -T +\frac{(1-\phi ^2)\alpha _1^2}{\sigma _\eta ^2} +\frac{1}{\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t^2\\&\quad +\frac{\rho }{\sigma _\eta (1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t\widetilde{R}_{t-1} -\left( a_{\sigma _\eta } +2\right) +\frac{b_{\sigma _\eta }}{\sigma _\eta ^2} \\ \nabla _{\tilde{\rho }}{\mathcal {L}}({\varSigma })&= (T-1)\rho -\frac{\rho }{\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{\alpha }_t^2 +\frac{1}{\sigma _\eta }\sum _{t=2}^{T}\widetilde{\alpha }_t\widetilde{R}_{t-1} \\&\quad +(a_\rho -1)(1-\rho ) -(b_\rho -1)(1+\rho ) \end{aligned}$$

To apply the leapfrog RMHMC algorithm, then we also need expressions for the metric tensor and its partial derivatives with respect to the transformed parameter. The metric tensor $\mathbf M (\sigma _\epsilon ,\sigma _\eta ,\rho )$ is a symmetric tridiagonal matrix, where the nonzero elements are calculated as follows:

$$\begin{aligned} \mathbf M (1,1)&= 2T +\frac{(T-1)\rho ^2}{1-\rho ^2} +\frac{2b_{\sigma _\epsilon }}{\sigma _\epsilon ^2}, \\ \mathbf M (2,2)&= 2T +\frac{(T-1)\rho ^2}{1-\rho ^2} +\frac{2b_{\sigma _\eta }}{\sigma _\eta ^2}, \\ \mathbf M (3,3)&= (T-1)(1+\rho ^2) +(a_\rho +b_\rho -2)(1-\rho ^2), \\ \mathbf M (1,2)&= \mathbf M (2,1) = -(T-1)\frac{\rho ^2}{1-\rho ^2},\\ \mathbf M (1,3)&= \mathbf M (3,1) = -(T-1)\rho ,\\ \mathbf M (2,3)&= \mathbf M (3,2) = -(T-1)\rho , \end{aligned}$$

and its partial derivatives follow as

$$\begin{aligned} \frac{\partial \mathbf M ({\varSigma })}{\partial \tilde{\sigma _\epsilon }}&= \left[ \begin{array}{ccc} -\dfrac{4b_{\sigma _\epsilon }}{\sigma _\epsilon ^2} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \end{array} \right] , ~~~ \frac{\partial \mathbf M ({\varSigma })}{\partial \tilde{\sigma _\eta }} = \left[ \begin{array}{ccc} 0 &{} 0 &{} 0 \\ 0 &{} -\dfrac{4b_{\sigma _\eta }}{\sigma _\eta ^2} &{} 0 \\ 0 &{} 0 &{} 0 \end{array} \right] , \\ \frac{\partial \mathbf M ({\varSigma })}{\partial \tilde{\rho }}&= -(T-1)\left[ \begin{array}{ccc} -\dfrac{2\rho }{1-\rho ^2} &{} \dfrac{2\rho }{1-\rho ^2} &{} 1-\rho ^2 \\ \dfrac{2\rho }{1-\rho ^2} &{} -\dfrac{2\rho }{1-\rho ^2} &{} 1-\rho ^2 \\ 1-\rho ^2 &{} 1-\rho ^2 &{} \dfrac{2(a_\rho +b_\rho -T-1)\rho (1-\rho ^2)}{T-1} \end{array} \right] . \end{aligned}$$

1.4 Updating parameters $\sigma _u^2$ and $c$

In the steps 4 and 5, we sample from the conditional posterior density of $\sigma _u^2$ and $c$ directly:

$$\begin{aligned} \sigma _u^2|c,\varvec{\alpha },\mathbf {Y}&\sim \mathcal {IG}\left( d_{\sigma _u}/2,D_{\sigma _u}/2\right) , \\ c|\sigma _u^2,\varvec{\alpha },\mathbf {Y}&\sim {\mathcal {N}}\left( M_c,V_c\right) , \end{aligned}$$

where

$$\begin{aligned}&d_{\sigma _u} = a_{\sigma _u}+T, ~~~ D_{\sigma _u} = b_{\sigma _u}+\sum \limits _{t=1}^{T} \left( Y_t-c-\alpha _t\right) ^2, \\&V_c = \left[ \frac{1}{v_c}+\frac{T}{\sigma _u^2}\right] ^{-1}, ~~~\hbox {and}~~~ M_c =V_c\left[ \frac{m_c}{v_c}+\frac{1}{\sigma _u^2}\sum \limits _{t=1}^{T} \left( Y_t-\alpha _t\right) \right] . \end{aligned}$$

Appendix 2: Estimation of RASVt model

Let us take the logarithm of the full conditional posterior density defined in Eq. (4):

$$\begin{aligned} \ln p(\theta _2,\mathbf h |\mathbf {R},\mathbf {Y})&= constant +\sum _{t=1}^{T}\left( -\frac{1}{2}h_t -\frac{1}{2}\ln (z_t) - \frac{1}{2}\widetilde{R}_t^2 \right) -T\ln (\sigma _\epsilon \sigma _u)\nonumber \\&\quad -\frac{1}{2\sigma _u^2}\sum _{t=1}^{T}\left( Y_t-\xi -h_t \right) ^2 -T\ln (\sigma _\eta ) + \frac{1}{2}\ln (1-\phi ^2) \nonumber \\&\quad -\frac{T-1}{2}\ln (1-\rho ^2)-\frac{1-\phi ^2}{2\sigma _\eta ^2}\widetilde{h}_1^2 -\frac{1}{2\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{h}_t^2\nonumber \\&\quad +\frac{\nu }{2}T\ln \left( \frac{\nu }{2}\right) -T\ln \varGamma \left( \frac{\nu }{2}\right) -\sum \nolimits _{t=1}^{T}\left[ \left( \frac{\nu }{2}+1\right) \ln (z_t)+\frac{\nu }{2} z_t^{-1}\right] \nonumber \\&\quad -\frac{1}{2}\ln (v_\xi ) -\frac{(\xi -m_\xi )^2}{2v_\xi } -\left( \frac{a_{\sigma _u}}{2}+1\right) \ln \left( \sigma _u^2\right) -\frac{b_{\sigma _u}}{2\sigma _u^2} \nonumber \\&\quad -\frac{1}{2}\ln (v_\mu ) -\frac{(\mu -m_\mu )^2}{2v_\mu } +(a_\phi -1)\ln (1+\phi ) +(b_\phi -1)\ln (1-\phi ) \nonumber \\&\quad -\left( \frac{a_{\sigma _\eta }}{2}+1\right) \ln \left( \sigma _\eta ^2\right) -\frac{b_{\sigma _\eta }}{2\sigma _\eta ^2} +(a_\rho -1)\ln (1+\rho )\nonumber \\&\quad +(b_\rho -1)\ln (1-\rho ) +(a_\nu -1)\ln (\nu ) -b_\nu \nu \end{aligned}$$

(6)

where

$$\begin{aligned} \widetilde{R}_t&= R_t\sigma _\epsilon ^{-1} e^{-\frac{1}{2}h_t} z_t^{-\frac{1}{2}},~~~\hbox {for}~t=1,\ldots ,T,~\hbox {and}\end{aligned}$$

(7)

$$\begin{aligned} \widetilde{h}_t&= \left\{ \begin{array}{l@{\quad }l} h_1 - \mu , &{} \text{ for } t=1 ;\\ h_t - \mu -\phi (h_{t-1}-\mu ) -\rho \sigma _\eta \widetilde{R}_{t-1}, &{} \text{ for } t=2,\ldots ,T.\end{array}\right. \end{aligned}$$

(8)

On the basis of the above density, the full conditional densities of $\mathbf h$, $\xi$, $\sigma _u$, $\phi$, and $\sigma _\eta$ are similar to $\varvec{\alpha }$, $c$, $\sigma _u$, $\phi$, and $\sigma _\eta$, respectively, for the RASV model. Thus, we follow the same sampling procedure described in the previous section with different specifications of $\widetilde{R}_t$ and $\widetilde{h}_t$ defined in Eqs. (7) and (8), respectively. To sample parameters $\mu$ and $\nu$ and latent variable $\mathbf z$, we describe these sampling procedures as follows.

1.1 Updating parameter $\mu$

On the basis of Eq. (6), parameter $\mu$ can be sampled directly from its conditional posterior density:

$$\begin{aligned} \mu |\phi ,\sigma _\eta ,\rho ,\mathbf h ,\mathbf z ,\mathbf R&\sim {\mathcal {N}}\left( M_{\mu },V_{\mu }\right) , \end{aligned}$$

where

$$\begin{aligned} V_{\mu }&= \left[ \frac{1}{v_{\mu }}+\frac{(1-\rho ^{2})(1-\phi ^{2})+(T-1)(1-\phi )^{2}}{\sigma _\eta ^2 (1-\rho ^2)}\right] ^{-1},\\ M_{\mu }&= V_{\mu }\left[ \frac{m_{\mu }}{v_{\mu }}+\frac{(1-\rho ^{2})(1-\phi ^2)h_1 +(1-\phi )\sum \limits _{t=2}^{T}(h_t-\phi h_{t-1}-\rho \sigma _\eta \widetilde{R}_{t-1})}{\sigma _\eta ^2 (1-\rho ^2)} \right] . \end{aligned}$$

1.2 Updating parameter $\nu$

The log posterior of $\nu$ and its gradient vector required in both leapfrog-HMC and -RMHMC algorithms are given by

$$\begin{aligned} {\mathcal {L}}(\nu )&= constant +\frac{\nu }{2}T\ln \left( \frac{\nu }{2}\right) -T\ln \varGamma \left( \frac{\nu }{2}\right) -\frac{\nu }{2}\sum \nolimits _{t=1}^{T}\left[ \ln (z_t)+z_t^{-1}\right] \\&\quad +(a_\nu -1)\ln (\nu ) -b_\nu \nu \end{aligned}$$

and

$$\begin{aligned} \nabla _\nu {\mathcal {L}}(\nu ) =\frac{1}{2}T\left[ \ln \left( \frac{\nu }{2}\right) +1\right] - \frac{1}{2}T\varPsi \left( \frac{\nu }{2}\right) - \frac{1}{2}\sum \limits _{t=1}^{T}\left[ \ln (z_t)+z_t^{-1}\right] +\frac{a_\nu -1}{\nu } - b_\nu . \end{aligned}$$

respectively, where $\varPsi (x)$ is a digamma function defined by $\varPsi (x)=d\ln \varGamma (x)/dx$. In the implementation of leapfrog-RMHMC algorithm, the metric tensor and its partial derivatives are calculated as follows:

$$\begin{aligned} \mathbf M (\nu ) = -\frac{T}{2\nu } + \frac{T}{4}\varPsi '\left( \frac{\nu }{2}\right) + \frac{a_\nu -1}{\nu ^2} \end{aligned}$$

and

$$\begin{aligned} \frac{\partial \mathbf M }{\partial \nu } = \frac{T}{2\nu ^2} + \frac{T}{8}\varPsi ''\left( \frac{\nu }{2}\right) -\frac{2(a_\nu -1)}{\nu ^3}, \end{aligned}$$

respectively, where $\varPsi '(x)$ is a trigamma function defined by $\varPsi '(x)=d\varPsi (x)/dx$ and $\varPsi ''(x)$ is a tetragamma function defined by $\varPsi ''(x)=d\varPsi '(x)/dx$.

1.3 Updating the latent variable z

The logarithm of the full conditional posterior density of the latent variable $\mathbf z$ is

$$\begin{aligned} {\mathcal {L}}(\mathbf z ) = constant +\sum _{t=1}^{T}\left[ -\frac{1}{2}\widetilde{R}_t^2 -\frac{\nu +3}{2}\ln (z_t) -\frac{\nu }{2z_t} \right] -\frac{1}{2\sigma _\eta ^2(1-\rho ^2)}\sum _{t=2}^{T}\widetilde{h}_t^2. \end{aligned}$$

From this conditional density, draws can be obtained by implementing the transformation $z_t=\exp (\hat{z}_t)$ to ensure constrained sampling for $z_t>0$. Gradient vector of the log posterior required in both leapfrog-HMC and -RMHMC algorithms is then evaluated as follows:

$$\begin{aligned} \nabla _{\hat{z}_t}{\mathcal {L}}(z_t) = - \widetilde{R}_t \frac{\partial \widetilde{R}_t}{\partial \hat{z}_t} - \frac{\nu +3}{2} + \frac{\nu }{2z_t} + \frac{\rho }{\sigma _\eta (1-\rho ^2)}\widetilde{h}_{t+1} \frac{\partial \widetilde{R}_t}{\partial \hat{z}_t} I_{t<T}, \end{aligned}$$

where $I$ is an indicator function and

$$\begin{aligned} \frac{\partial \widetilde{R}_t}{\partial \hat{z}_t} = -\dfrac{1}{2}R_t \sigma _\epsilon ^{-1}e^{-\frac{1}{2}h_t}z_t^{-\frac{1}{2}}. \end{aligned}$$

Furthermore, the metric tensor $\mathbf M (\mathbf z )$ required to apply the leapfrog-RMHMC algorithm is a diagonal matrix whose diagonal entries are given by

$$\begin{aligned} \mathbf M (t,t) = \dfrac{\nu +1}{2} + \dfrac{\rho ^2}{4(1-\rho ^2)}I_{t<T}. \end{aligned}$$

Since the above metric tensor is fixed, that is not a function of the latent variable z, the associated partial derivatives with respect to the transformed latent variable are zero. This implies that both metric tensor and its inverse are not recalculated in the generalized leapfrog integration scheme.