An application of an objective prior for the number of degrees of freedom of a multivariate t distribution
F. Javier Rubio and Cristiano Villa
14 July, 2017
A financial application of an objective prior for the number of degrees of freedom of a multivariate t distribution
The \(d\)-variate \(t\) probability density function with \(\nu>0\) degrees of freedom (see and for an extensive review of this model) is given by \[ f_d({\bf x}\vert {\mu},{\Sigma},\nu) = \dfrac{\Gamma\left(\dfrac{\nu+d}{2}\right)}{\Gamma\left(\dfrac{\nu}{2}\right)\sqrt{(\pi\nu)^d\vert{\Sigma}\vert}} \left(1+\dfrac{({\bf x}-{\mu})^{T}{\Sigma}^{-1}({\bf x}-{\mu})}{\nu}\right)^{-\frac{\nu+d}{2}}, \] [1] derived a discrete objective prior for the number of degrees of freedom of the multivariate \(t\) distribution. This prior is truncated at \(\nu_{\max}\), pre-specified by the user, typically \(\nu_{\max}=30\) and discrete, see [3] for more details.
Here, we present a financial application of the multivariate \(t\) distribution, using the prior for the degrees of freedom parameter proposed in [1], where we model jointly the daily log-returns for the Swiss Market Index (SMI) and Swiss reinsurer (Swiss.Re). The data are available from the R package ‘ghyp’ and contain \(n = 1769\) observations corresponding to the period January 2000 to January 2007.
References
Kullback Leibler divergence between two multivariate t distributions
An objective prior for the number of degrees of freedom of a multivariate t distribution
##############################################################################################
# Data from : Statistics and Data Analysis for Financial Engineering pp. 180
##############################################################################################
rm(list=ls())
library(mnormt)
library(mvtnorm)
library(mcmc)
library(Ecdat)## Loading required package: Ecfun## Warning: package 'Ecfun' was built under R version 3.3.3##
## Attaching package: 'Ecfun'## The following object is masked from 'package:base':
##
## sign##
## Attaching package: 'Ecdat'## The following object is masked from 'package:datasets':
##
## Orangelibrary(knitr)## Warning: package 'knitr' was built under R version 3.3.3data(CRSPday, package = "Ecdat")
Y = CRSPday[ , c(5, 7)]
##############################################################################################
# Optimisation step as in the book
##############################################################################################
loglik = function(par){
mu = par[1:2]
A = matrix(c(par[3], par[4], 0, par[5]), nrow = 2, byrow = T)
scale_matrix = t(A) %*% A
df = par[6]
-sum(log(dmt(Y, mean = mu, S = scale_matrix, df = df)))
}
A = chol(cov(Y))
start = as.vector(c(apply(Y, 2, mean),A[1, 1], A[1, 2], A[2, 2], 4))
fit_mvt = optim(start, loglik, method = "L-BFGS-B",lower = c(-0.02, -0.02, -0.1, -0.1, -0.1, 2),upper = c(0.02, 0.02, 0.1, 0.1, 0.1, 15), hessian = T)
##############################################################################################
# Optimisation step
##############################################################################################
llt = function(par){
if(par[3]>0 & par[5]>0 & par[6]>0){
mu = par[1:2]
A = matrix(c(par[3], par[4], 0, par[5]), nrow = 2, byrow = T)
Sigma = t(A) %*% A
var = dmvt(Y, delta = mu, sigma = Sigma, df = par[6], log = TRUE,type = "shifted")
return( -sum(var) )
}
else return(Inf)
}
loglik(start)## [1] -15839.85llt(start)## [1] -15839.85OPT = optim(start,llt,control=list(maxit=1000))
# The optimisation implemented here is slightly better
OPT$val## [1] -16106.81fit_mvt$val## [1] -16105.65OPT$par## [1] 0.0004990631 0.0008414805 0.0125729426 0.0026557420 0.0049487891
## [6] 4.1890062126fit_mvt$par## [1] 0.0003789236 0.0008317070 0.0126907363 0.0026859493 0.0051011198
## [6] 4.2618395311##############################################################################################
# log prior
##############################################################################################
# Kullback Leibler in two dimensions
rm(d)## Warning in rm(d): object 'd' not foundd=2 # dimension
# Normalising constant
K <- function(d,nu) (gamma(0.5*(nu+d))/( gamma(0.5*nu)*sqrt((pi*nu)^d) ))
# Kullback Liebler divergence
DKL <- function(nu,nup){
val1 <- log(K(d,nu)/K(d,nup)) -0.5*(nu+d)*(digamma(0.5*(nu+d))-digamma(0.5*nu))
tempf <- Vectorize(function(t) (1+t/nu)^(-0.5*(nu+d))*t^(0.5*d-1)*log(1+t/nup))
int<- integrate(tempf,0,Inf,rel.tol = 1e-9)$value
val2 <- (K(d,nu)*pi^(0.5*d)*0.5*(nup+d)/gamma(0.5*d))*int
return(val1+val2)
}
numax = 30
T2 <- matrix(0,ncol=3,nrow=numax)
T2[1,] <- c(1,0,DKL(1,2))
for(i in 2:numax) T2[i,] = c(i,DKL(i,i-1),DKL(i,i+1))
colnames(T2) <- c("nu","DKL(nu,nu-1)","DKL(nu,nu+1)")
rownames(T2) <- 1:numax
print(kable(T2,digits=12))##
##
## nu DKL(nu,nu-1) DKL(nu,nu+1)
## --- ------------- -------------
## 1 0.000000e+00 1.415927e-01
## 2 7.944154e-02 2.732554e-02
## 3 1.956127e-02 9.139055e-03
## 4 7.282898e-03 3.961126e-03
## 5 3.352862e-03 2.005447e-03
## 6 1.763603e-03 1.127459e-03
## 7 1.017649e-03 6.838292e-04
## 8 6.288552e-04 4.393997e-04
## 9 4.097277e-04 2.954756e-04
## 10 2.784711e-04 2.061445e-04
## 11 1.959103e-04 1.482667e-04
## 12 1.418505e-04 1.094059e-04
## 13 1.052417e-04 8.251520e-05
## 14 7.973102e-05 6.342075e-05
## 15 6.151054e-05 4.955533e-05
## 16 4.821469e-05 3.928773e-05
## 17 3.832778e-05 3.155168e-05
## 18 3.085192e-05 2.563259e-05
## 19 2.511425e-05 2.104082e-05
## 20 2.065128e-05 1.743421e-05
## 21 1.713759e-05 1.456933e-05
## 22 1.434074e-05 1.227021e-05
## 23 1.209212e-05 1.040777e-05
## 24 1.026762e-05 8.886050e-06
## 25 8.774727e-06 7.632841e-06
## 26 7.543654e-06 6.593187e-06
## 27 6.521166e-06 5.724832e-06
## 28 5.666240e-06 4.994967e-06
## 29 4.946970e-06 4.377892e-06
## 30 4.338319e-06 3.853307e-06prior = Vectorize(function(nu){
var =ifelse(nu<(numax-1), exp( T2[nu,3]),exp( T2[nu,2]))
return(var-1)
})
vec = 1:numax
pr1 = prior(vec)
npr1 = pr1/sum(pr1)
lprior = Vectorize(function(nu) log(npr1[nu]))
# Objective prior for nu
plot(npr1,ylim=c(0,0.8),xlab=~nu,ylab="Prior probability",main="Objective prior",cex.axis=1.5,cex.lab=1.5,pch=19,cex=1.5)
##############################################################################################
# log posterior
##############################################################################################
lpt = function(par){
if(par[3]>0 & par[5]>0 & par[6]>0 & par[6]<31 & par[6]%%1==0){
mu = par[1:2]
A = matrix(c(par[3], par[4], 0, par[5]), nrow = 2, byrow = T)
Sigma = t(A) %*% A
var = dmvt(Y, delta = mu, sigma = Sigma, df = par[6], log = TRUE,type = "shifted")
return( sum(var) + lprior(par[6]) - 1.5*log(det(Sigma)) )
}
else return(-Inf)
}
# Function to sample from the discrete parameter nu
post.samp.nu =function(par) {
lpr = probr = vector()
for(k in 1:30) lpr[k] = lpt(c(par[1:5],k))
for(k in 1:30) probr[k] = 1/sum(exp(lpr-lpr[k]))
return(sample(1:30, 1, prob = probr))
}
# number of iterations
NS= 55000
# Sampler (Metropolis step for the continuous parameters)
# Approximately 10 minutes running time (status bar available)
samp = matrix(0,nrow=NS,ncol=6)
samp[1,] = start
# Initiate the status bar
#sb <- txtProgressBar(min = 0, max = NS, style = 3)
for(i in 2:NS){
samp[i,6] = post.samp.nu(samp[i-1,1:5])
tempf = function(par) lpt(c(par,samp[i,6]))
samp[i,1:5] = metrop(tempf, samp[i-1,1:5], 1, scale = c(0.0003,0.0001,0.0002,0.0001,0.0001))$final
# Update the status var
# setTxtProgressBar(sb, i)
}
# Some histograms and traceplots
burn = 5000
thin = 50
index = seq(burn,NS,thin)
mu1p = samp[index,1]
mu2p = samp[index,2]
sigma1p = samp[index,3]
sigma2p = samp[index,4]
sigma3p = samp[index,5]
nup = samp[index,6]
hist(mu1p)
hist(mu2p)
hist(sigma1p)
hist(sigma2p)
hist(sigma3p)
hist(nup)
plot(mu1p,type="l")
plot(mu2p,type="l")
plot(sigma1p,type="l")
plot(sigma2p,type="l")
plot(sigma3p,type="l")
plot(nup,type="l")
# Comparing posterior means with MLE
colMeans(samp[index,])## [1] 0.0004301167 0.0008573191 0.0124271193 0.0026263064 0.0049170305
## [6] 4.0079920080OPT$par## [1] 0.0004990631 0.0008414805 0.0125729426 0.0026557420 0.0049487891
## [6] 4.1890062126###############################################################################
# Posterior density plots using the posterior mean estimators
###############################################################################
postmean = colMeans(samp[index,])
post.mean.den = function(x){
mu = postmean[1:2]
A = matrix(c(postmean[3], postmean[4], 0, postmean[5]), nrow = 2, byrow = T)
Sigma = t(A) %*% A
var = dmvt(x, delta = mu, sigma = Sigma, df = postmean[6], type = "shifted",log=FALSE)
return(var)
}
library(fMultivar)## Loading required package: timeDate## Loading required package: timeSeries## Loading required package: fBasics## ## Rmetrics Package fBasics## Analysing Markets and calculating Basic Statistics## Copyright (C) 2005-2014 Rmetrics Association Zurich## Educational Software for Financial Engineering and Computational Science## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.## https://www.rmetrics.org --- Mail to: info@rmetrics.org## ## Rmetrics Package fMultivar## Analysing and Modeling Multivariate Financial Return Distributions## Copyright (C) 2005-2014 Rmetrics Association Zurich## Educational Software for Financial Engineering and Computational Science## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.## https://www.rmetrics.org --- Mail to: info@rmetrics.org# Bivariate fitted posterior density
x <- (-100:100)/1000
X <- grid2d(x)
z <- apply(cbind(X$x,X$y),1,post.mean.den)
ZD <- list(x = x, y = x, z = matrix(z, ncol = length(x)))
# Perspective Density Plot:
persp(ZD, theta = 60, phi = 20, col = "steelblue")
# Contour Density Plot:
plot(as.vector(Y[,1]),as.vector(Y[,2]),xlim=c(-0.15,0.15),ylim=c(-0.075,0.05), main="Fitted Posterior Density Contours",xlab="IBM",ylab="CRSP",cex.axis=1.5,cex.lab=1.5)
par(new=TRUE)
contour(ZD,levels = c(0.55,1,2,4,8,16,32,64,128,256,512,1024,2048),xlim=c(-0.15,0.15),ylim=c(-0.075,0.05),col="red",cex.axis=1.5,cex.lab=1.5,lwd=2)