SV with Gaussian or Student’s \(t\) errors

Hamiltonian MC via STAN

1 Downloading stock returns data

#install.packages("rstan")
#install.packages("astsa")
library("rstan")
library("astsa")
options(mc.cores=4)

# Daily closing prices of Germany DAX
price = EuStockMarkets[,4]
y     = diff(log(price))
y     = y - mean(y)
n     = length(y)

par(mfrow=c(1,1))
ts.plot(y,ylab="returns",main="Germany DAX")

2 SV-AR(1) model

For \(t=1,\ldots,n\), \[\begin{eqnarray*} y_t &=& \exp\{h_t/2\}\varepsilon_t\\ h_t &=& \mu+\phi (h_{t-1} - \mu)+\sigma u_t,\\ h_1 &\sim& N(\mu,\sigma^2/(1-\phi^2)). \end{eqnarray*}\] where \(\varepsilon_t\) are iid \(t_\nu(0,1)\) and \(u_t\) are iid \(N(0,1)\). The independent priors used are: \(\mu \sim N(0,10^2)\), \(\phi \sim N_{(-1,1)}(0.8,0.3^2)\), \(\sigma^2 \sim Gamma(0.5,0.5)\) and \(\nu \sim Gamma(1,0.05)\).

3 Stan model

3.1 Gaussian errors

# Filename: sv-g.stan
data {
  int<lower=0> n;   // # time points (equally spaced)
  vector[n] y;      // mean corrected return at time t
}

parameters {
  real mu;                     // mean log volatility
  real<lower=-1,upper=1> phi;  // persistence of volatility
  real<lower=0> sigma;         // white noise shock scale
  real<lower=2> nu;            // degrees of freedom
  vector[n] h;                 // log volatility at time t
}

model {
  phi ~ normal(0.8,0.3);
  sigma ~ gamma(0.5,0.5);
  mu ~ normal(0,10);
  nu ~ gamma(1,0.05);
  h[1] ~ normal(mu,sigma/sqrt(1-phi*phi));
  for (t in 2:n)
    h[t] ~ normal(mu+phi*(h[t-1]-mu),sigma);
  y ~ student_t(nu, 0.0, exp(h/2));
}

3.2 Student’s \(t\) errors

# Filename: sv-t.stan
data {
  int<lower=0> n;   // # time points (equally spaced)
  vector[n] y;      // mean corrected return at time t
}

parameters {
  real mu;                     // mean log volatility
  real<lower=-1,upper=1> phi;  // persistence of volatility
  real<lower=0> sigma;         // white noise shock scale
  vector[n] h;                 // log volatility at time t
}

model {
  phi ~ normal(0.8,0.3);
  sigma ~ gamma(0.5,0.5);
  mu ~ normal(0,10);
  h[1] ~ normal(mu,sigma/sqrt(1-phi*phi));
  for (t in 2:n)
    h[t] ~ normal(mu+phi*(h[t-1]-mu),sigma);
  y ~ normal(0.0,exp(h/2));
}

4 Running stan code

M  = 500

model.n = stan_model("sv-g.stan")

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

fit.g   = sampling(model.n,list(n=n,y=y),iter=M,chains=1)

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000226 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 6.575 seconds (Warm-up)
## Chain 1:                7.361 seconds (Sampling)
## Chain 1:                13.936 seconds (Total)
## Chain 1:

model.t = stan_model("sv-t.stan")

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

fit.t   = sampling(model.t,list(n=n,y=y),iter=M,chains=1)

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000303 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.03 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1: Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1: Iteration: 500 / 500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 7.734 seconds (Warm-up)
## Chain 1:                6.858 seconds (Sampling)
## Chain 1:                14.592 seconds (Total)
## Chain 1:

params.g = extract(fit.g)
params.t = extract(fit.t)

qsd.g = t(apply(exp(params.g$h/2),2,quantile,c(0.025,0.5,0.975)))
qsd.t = t(apply(exp(params.t$h/2),2,quantile,c(0.025,0.5,0.975)))

5 Marginal posteriors of \((\mu,\phi,\sigma,\nu)\)

par(mfrow=c(2,4))
ts.plot(params.g$mu,xlab="Draws",main=expression(mu),ylab="",ylim=range(params.g$mu,params.t$mu))
lines(params.t$mu,col=2)
ts.plot(params.g$phi,xlab="Draws",main=expression(phi),ylab="",ylim=range(params.g$phi,params.t$phi))
lines(params.t$phi,col=2)
ts.plot(params.g$sigma,xlab="Draws",main=expression(sigma),ylab="",ylim=range(params.g$sigma,params.t$sigma))
lines(params.t$sigma,col=2)
ts.plot(params.t$nu,xlab="Iterations",main=expression(nu),ylab="")
boxplot(params.g$mu,params.t$mu,outline=FALSE,names=c("Gaussian","Student's t"),main=expression(mu))
boxplot(params.g$phi,params.t$phi,outline=FALSE,names=c("Gaussian","Student's t"),main=expression(phi))
boxplot(params.g$sigma,params.t$sigma,outline=FALSE,names=c("Gaussian","Student's t"),main=expression(sigma))
hist(params.t$nu,prob=TRUE,main=expression(nu),xlab="")

5.1 Posterior quantiles for \(\mu\)

tab = rbind(quantile(params.g$mu,c(0.01,0.05,0.5,0.95,0.99)),
quantile(params.t$mu,c(0.01,0.05,0.5,0.95,0.99)))
rownames(tab)=c("Gaussian","Student's t")
tab

##                    1%        5%       50%       95%       99%
## Gaussian    -10.07859 -10.01688 -9.827687 -9.642983 -9.538947
## Student's t -10.16613 -10.09410 -9.930063 -9.748884 -9.667854

5.2 Posterior quantiles for \(\phi\)

tab = rbind(quantile(params.g$phi,c(0.01,0.05,0.5,0.95,0.99)),
quantile(params.t$phi,c(0.01,0.05,0.5,0.95,0.99)))
rownames(tab)=c("Gaussian","Student's t")
tab

##                    1%        5%       50%       95%       99%
## Gaussian    0.9437874 0.9529665 0.9664373 0.9775090 0.9804357
## Student's t 0.9358860 0.9401383 0.9661100 0.9804359 0.9844241

5.3 Posterior quantiles for \(\sigma\)

tab = rbind(quantile(params.g$sigma,c(0.01,0.05,0.5,0.95,0.99)),
quantile(params.t$sigma,c(0.01,0.05,0.5,0.95,0.99)))
rownames(tab)=c("Gaussian","Student's t")
tab

##                    1%        5%       50%       95%       99%
## Gaussian    0.1306978 0.1364088 0.1518336 0.1813667 0.1898979
## Student's t 0.1149745 0.1183771 0.1528463 0.1946473 0.2083977

5.4 Posterior quantiles for \(\nu\)

quantile(params.t$nu,c(0.01,0.05,0.5,0.95,0.99))

##       1%       5%      50%      95%      99% 
## 10.00972 11.68812 20.81765 46.45052 66.84300

6 Posterior distribution of time-varying standard deviations

par(mfrow=c(1,1))
ts.plot(qsd.g[,2],ylim=range(qsd.g,qsd.t),ylab="Standard deviation")
lines(qsd.g[,1])
lines(qsd.g[,3])
lines(qsd.t[,1],col=2)
lines(qsd.t[,2],col=2)
lines(qsd.t[,3],col=2)
legend("top",legend=c("Gaussian errors","Student's t errors"),col=1:2,lwd=2,bty="n",lty=1)