Petrobras returns

Gaussian, Student’s \(t\) or SV?

1 Petrobras data

#install.packages("tidyquant")
library("tidyquant")

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## ── Attaching core tidyquant packages ──────────────────────── tidyquant 1.0.9 ──
## ✔ PerformanceAnalytics 2.0.4      ✔ TTR                  0.24.4
## ✔ quantmod             0.4.26     ✔ xts                  0.13.2
## ── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
## ✖ zoo::as.Date()                 masks base::as.Date()
## ✖ zoo::as.Date.numeric()         masks base::as.Date.numeric()
## ✖ PerformanceAnalytics::legend() masks graphics::legend()
## ✖ quantmod::summary()            masks base::summary()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

getSymbols('PBR',src='yahoo',return.class='ts',from="2019-01-01")

## [1] "PBR"

n = nrow(PBR)

ret = log(PBR[2:n,6]/PBR[1:(n-1),6])

y   = ret-mean(ret)
n   = n - 1

par(mfrow=c(1,1))
ts.plot(PBR[,4],main="",ylab="Price")

par(mfrow=c(1,2))
ts.plot(ret,main="",ylab="Log-return")
hist(ret,breaks=seq(-max(abs(ret)),max(abs(ret)),length=50),prob=TRUE,xlab="Log-return",main="")

2 Model 1: iid Gaussian

y_t iid N(0,sig2)

sig.hat.1 = sqrt(mean(y^2))
sig.hat.1

## [1] 0.03279088

3 Gaussian vs Student’s \(t\)

xxx = seq(-10,10,length=1000)
par(mfrow=c(1,1))
plot(xxx,dnorm(xxx),type="l",xlab="x",ylab="Density")
lines(xxx,dt(xxx,df=1),col=2)
lines(xxx,dt(xxx,df=2),col=3)
lines(xxx,dt(xxx,df=30),col=4)
legend("topleft",legend=c("df=1","df=2","df=30","Gaussian"),col=c(2:4,1),lty=1,lwd=2)

par(mfrow=c(1,1))
plot(xxx,log(dnorm(xxx)),type="l",xlab="x",ylab="Log density")
lines(xxx,log(dt(xxx,df=1)),col=2)
lines(xxx,log(dt(xxx,df=2)),col=3)
lines(xxx,log(dt(xxx,df=30)),col=4)

4 Model 2: iid Student’s \(t\)

sigs = seq(0.001,0.15,length=1000)
nus  = c(1:30,seq(35,100,by=5),seq(200,1000,by=100))
nnu  = length(nus)
like = rep(0,1000)
sig.mle = rep(0,nnu)
likes = rep(0,nnu)
for (j in 1:nnu){
  for (i in 1:1000)
    like[i] = sum(dt(y/sigs[i],df=nus[j],log=TRUE))-n*log(sigs[i])
  sig.mle[j] = sigs[like==max(like)]
  likes[j] = like[like==max(like)]
}

sig.hat.2 = sig.mle[likes==max(likes)]
nu.hat = nus[likes==max(likes)]

c(sig.hat.1,sig.hat.2,nu.hat)

## [1] 0.03279088 0.01934535 3.00000000

par(mfrow=c(1,1))
plot(nus,sig.mle,ylim=range(sig.mle,sig.hat.1))
abline(h=sig.hat.1)

par(mfrow=c(1,1))
plot(likes,sig.mle,col=0)
text(likes,sig.mle,nus,cex=0.5)

par(mfrow=c(1,1))
xxx = seq(-max(abs(ret)),max(abs(ret)),length=1000)
xbreaks = seq(-max(abs(ret)),max(abs(ret)),length=50)
hist(ret,breaks=xbreaks,prob=TRUE,xlab="Log-return",main="")
lines(xxx,dnorm(xxx,0,sig.hat.1),col=2,lwd=3)
lines(xxx,dt(xxx/sig.hat.2,df=nu.hat)/sig.hat.2,col=4,lwd=2)

par(mfrow=c(1,1))
legends  = c("Data","Gaussian model","Student's t model")
alphas   = seq(0.01,0.99,by=0.02)
qdata    = quantile(ret,alphas)
qnormal  = qnorm(alphas,0,sig.hat.1)
qstudent = sig.hat.2*qt(alphas,df=nu.hat)
xrange   = range(qdata,qnormal,qstudent)
plot(qdata,alphas,type="l",xlim=xrange,ylab="Tail probabilities",xlab="Quantiles",lwd=2)
lines(qnormal,alphas,type="l",col=2,lwd=2)
lines(qstudent,alphas,type="l",col=4,lwd=2)
legend("topleft",legend=legends,col=c(1,2,4),lwd=3,lty=1)

5 Stochastic volatility models

Here we assume that \[\begin{eqnarray*} y_t &=& e^{h_t/2} \varepsilon_t\\ h_{t+1} &=& \mu + \phi (h_t-\mu) + \sigma u_t, \end{eqnarray*}\] with \(cor(\varepsilon_t,u_t)=\rho\).

5.1 SV, SVt, SVl, SVtl

Four models contemplated by Gregor Kastner’s R package stochvol:

• SV with Gaussian errors (SV) \[ \varepsilon_t \sim N(0,1),\ u_t \sim N(0,1), \ \rho=0 \]

• SV with Student’s \(t\) errors (SVt) \[ \varepsilon_t \sim t_\nu(0,1),\ u_t \sim N(0,1), \ \rho=0, \ \nu \ \mbox{unknown} \]

• SV with Gaussian errors and leverage effect (SVl) \[ \varepsilon_t \sim N(0,1),\ u_t \sim N(0,1), \ \rho \ \mbox{unknown}, \]

• SV with Student’s \(t\) errors and leverage effect (SVtl) \[ \varepsilon_t \sim N(0,1),\ u_t \sim N(0,1), \ \rho \ \mbox{unknown}, \ \nu \ \mbox{unknown} \]

5.2 Prior distributions

The package stochvol has the following default prior specifications: \[\begin{eqnarray*} \mu &\sim& N(0,100^2)\\ (\phi+1)/2 &\sim& Beta(5,1.5)\\ \sigma^2 &\sim& Gamma(1/2,1/2)\\ h_1 &\sim& N(\mu,\sigma^2/(1-\phi^2))\\ (\rho+1)/2 &\sim& Beta(4,4)\\ \nu &\sim& Exp(0.1) \end{eqnarray*}\]

5.3 References

• Source: Modeling Univariate and Multivariate Stochastic Volatility in R with stochvol and factorstochvol - https://cran.r-project.org/web/packages/factorstochvol/vignettes/paper.pdf

• Kastner (2016) Dealing with Stochastic Volatility in Time Series Using the R Package stochvol. Journal of Statistical Software, 69(5), 1-30. doi:10.18637/jss.v069.i05.

• Kastner(2019) Sparse Bayesian Time-Varying Covariance Estimation in Many Dimensions. Journal of Econometrics, 210(1), 98-115. doi:10.1016/j.jeconom.2018.11.007.

• Kastner and Fruhwirth-Schnatter (2014) Ancillarity-Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Estimation of Stochastic Volatility Models. Computational Statistics and Data Analysis, 76, 408–423. doi:10.1016/j.csda.2013.01.002.

• Kastner, Fruhwirth-Schnatter and Lopes (2017) Efficient Bayesian Inference for Multivariate Factor Stochastic Volatility Models. Journal of Computational and Graphical Statistics, 26(4), 905–917. doi:10.1080/10618600.2017.1322091.

library("stochvol")
M        = 10000
burn     = 10000
fit.sv   = svsample(ret,draws=M,burnin=burn,quiet=TRUE)

## Argument 'y' (data vector) contains values very close to zero. I am applying an offset constant of size 3.28023364663639e-06 to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.

fit.svt  = svtsample(ret,draws=M,burnin=burn,quiet=TRUE)

## Argument 'y' (data vector) contains values very close to zero. I am applying an offset constant of size 3.28023364663639e-06 to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.

fit.svl  = svlsample(ret,draws=M,burnin=burn,quiet=TRUE)

## Argument 'y' (data vector) contains values very close to zero. I am applying an offset constant of size 3.28023364663639e-06 to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.

fit.svtl = svtlsample(ret,draws=M,burnin=burn,quiet=TRUE)

## Argument 'y' (data vector) contains values very close to zero. I am applying an offset constant of size 3.28023364663639e-06 to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.

qlogvols = cbind(apply(fit.sv$latent[[1]],2,median),apply(fit.svt$latent[[1]],2,median),
apply(fit.svl$latent[[1]],2,median),apply(fit.svtl$latent[[1]],2,median))

qstd.sv   = t(apply(exp(fit.sv$latent[[1]]/2),2,quantile,c(0.05,0.5,0.95)))
qstd.svt  = t(apply(exp(fit.svt$latent[[1]]/2),2,quantile,c(0.05,0.5,0.95)))
qstd.svl  = t(apply(exp(fit.svl$latent[[1]]/2),2,quantile,c(0.05,0.5,0.95)))
qstd.svtl = t(apply(exp(fit.svtl$latent[[1]]/2),2,quantile,c(0.05,0.5,0.95)))

par(mfrow=c(2,3))
pars = cbind(exp(fit.sv$para[[1]][,1]/2),exp(fit.svt$para[[1]][,1]/2),
             exp(fit.svl$para[[1]][,1]/2),exp(fit.svtl$para[[1]][,1]/2))
plot(density(pars[,1]),xlab=expression(mu),main="",xlim=c(0,0.05))
for (i in 2:4) lines(density(pars[,i]),col=i)
legend("topright",legend=c("SV","SVt","SVl","SVtl"),col=1:4,lty=1)

pars = cbind(fit.sv$para[[1]][,2],fit.svt$para[[1]][,2],fit.svl$para[[1]][,2],fit.svtl$para[[1]][,2])
plot(density(pars[,1]),xlab=expression(phi),main="",xlim=range(pars),ylim=c(0,40))
for (i in 2:4) lines(density(pars[,i]),col=i)

pars = cbind(fit.sv$para[[1]][,3],fit.svt$para[[1]][,3],fit.svl$para[[1]][,3],fit.svtl$para[[1]][,3])
plot(density(pars[,1]),xlab=expression(sigma),main="",xlim=range(pars),ylim=c(0,20))
for (i in 2:4) lines(density(pars[,i]),col=i)

pars = cbind(fit.svt$para[[1]][,4],fit.svtl$para[[1]][,4])
plot(density(pars[,1]),xlab=expression(nu),main="",xlim=range(pars),col=2,ylim=c(0,0.5))
lines(density(pars[,2]),col=4)

pars = cbind(fit.svl$para[[1]][,5],fit.svtl$para[[1]][,5])
plot(density(pars[,1]),xlab=expression(rho),main="",xlim=range(pars),col=3)
lines(density(pars[,2]),col=4)

par(mfrow=c(1,1))
plot(qlogvols[,1],type="l",ylim=range(qlogvols),xlab="Time",ylab="Log volatilities")
lines(qlogvols[,2],col=2)
lines(qlogvols[,3],col=3)
lines(qlogvols[,4],col=4)
legend("topleft",legend=c("SV","SVt","SVl","SVtl"),col=1:4,lty=1)

yrange = range(0,qstd.sv[,2],qstd.svt[,2],qstd.svl[,2],qstd.svtl[,2])
par(mfrow=c(1,1))
plot(qstd.sv[,2],ylab="Standard deviations",ylim=yrange,type="l",xlab="Time")
lines(qstd.svt[,2],col=2)
lines(qstd.svl[,2],col=3)
lines(qstd.svtl[,2],col=4)
legend("topleft",legend=c("SV","SVt","SVl","SVtl"),col=1:4,lty=1)

par(mfrow=c(2,2))
yrange = range(qstd.sv,qstd.svt,qstd.svl,qstd.svtl)
ts.plot(qstd.sv,ylab="Standard deviations",col=c(3,2,3),ylim=yrange)
lines(0.1*abs(ret),type="h");title("SV")
ts.plot(qstd.svt,ylab="Standard deviations",col=c(3,2,3),ylim=yrange)
lines(0.1*abs(ret),type="h");title("SVt")
ts.plot(qstd.svl,ylab="Standard deviations",col=c(3,2,3),ylim=yrange)
lines(0.1*abs(ret),type="h");title("SVl")
ts.plot(qstd.svtl,ylab="Standard deviations",col=c(3,2,3),ylim=yrange)
lines(0.1*abs(ret),type="h");title("SVtl")

par(mfrow=c(2,2))
for (i in 1:4){
  ind = ((i-1)*n/4+1):(i*n/4)
  yrange = range(0,qstd.sv[ind,2],qstd.svt[ind,2],qstd.svl[ind,2],qstd.svtl[ind,2])
  plot(ind,qstd.sv[ind,2],ylab="Standard deviations",ylim=yrange,type="l",xlab="Time")
  lines(ind,qstd.svt[ind,2],col=2)
  lines(ind,qstd.svl[ind,2],col=3)
  lines(ind,qstd.svtl[ind,2],col=4)
  lines(ind,0.3*abs(ret[ind]),type="h")
}
legend("topleft",legend=c("SV","SVt","SVl","SVtl"),col=1:4,lty=1)