Abstract
This vignette is a tutorial for a R packageinvgamstochvol
.
The main function lik_clo
computes the log likelihood for
an inverse gamma stochastic volatility model using a closed form
expression of the likelihood. The closed form expression is obtained for
the log likelihood of a stationary inverse gamma stochastic volatility
model by marginalising out the volatilities. This allows the user to
obtain the maximum likelihood estimator for this non linear non Gaussian
state space models. In addition, we can obtain the smoothed estimates of
the volatility using draws from the exact posterior distribution of the
inverse volatility by calling the function DrawK0
. Lastly
one can evaluate the 2F1 hypergeometric function using
ourgeo
. The computation of this closed form expression
details are given in Leon-Gonzalez, R., & Majoni, B. (2023). Exact
Likelihood for Inverse Gamma Stochastic Volatility Models (No. 23-11).
We first provide an overview of the package, and then a quick tutorial.
The closed form expression is obtained for the log likelihood of a stationary inverse gamma stochastic volatility model by marginalising out the volatilities. This allows the user to obtain the maximum likelihood estimator for this non linear non Gaussian state space model. Further, the package can provide the smoothed estimates of the volatility by averaging draws from the exact posterior distribution of the inverse volatilities.
#install.packages("invgamstochvol")
library(invgamstochvol)
The data set that we use for this example has 150 observations. Ydep are the observed data, rho represents the parameter for the persistence of the volatility, p is the number of lags and Xdep are the regressors.
##simulate data
n=150
dat<-data.frame(Ydep=runif(n,0.3,1.4))
Ydep <- as.matrix(dat, -1,ncol=ncol(dat))
littlerho=0.95
r0=1
rho=diag(r0)*littlerho
p=4
n=4.1
T=nrow(Ydep)
Xdep <- Ydep[p:(T-1),]
if (p>1){
for(lagi in 2:p){
Xdep <- cbind(Xdep, Ydep[(p-lagi+1):(T-lagi),])
}
}
T=nrow(Ydep)
Ydep <- as.matrix(Ydep[(p+1):T,])
T=nrow(Ydep)
unos <- rep(1,T)
Xdep <- cbind(unos, Xdep)
The matrix of residuals from OLS can be obtained as follows.
## obtain residuals
bOLS <- solve(t(Xdep) %*% Xdep) %*% t(Xdep) %*% Ydep
Res= Ydep- Xdep %*% bOLS
Res=Res[1:T,1]
b2=solve(t(Res) %*% Res/T) %*% (1-rho %*% rho)/(n-2)
Res=as.matrix(Res,ncol=1)
The function lik_clo returns a list of 7 items. List item number 1, is the sum of the log likelihood, while the rest are constants that are useful to obtain the smoothed estimates of the volatility.
## obtain the log likelihood
LL1=lik_clo(Res,b2,n,rho)
LL1[1]
## [[1]]
## [1] -40.16983
To obtain likelihood, the same approach as highlighted in the example using simulated data above applies. After obtaining the likelihood, we show how the smoothed estimates of volatility can be obtained.
##Example using US data
data1 <- US_Inf_Data
Ydep <- as.matrix(data1)
littlerho=0.95
r0=1
rho=diag(r0)*littlerho
p=4
n=4.1
T=nrow(Ydep)
Xdep <- Ydep[p:(T-1),]
if (p>1){
for(lagi in 2:p){
Xdep <- cbind(Xdep, Ydep[(p-lagi+1):(T-lagi),])
}
}
T=nrow(Ydep)
Ydep <- as.matrix(Ydep[(p+1):T,])
T=nrow(Ydep)
unos <- rep(1,T)
Xdep <- cbind(unos, Xdep)
## obtain residuals
bOLS <- solve(t(Xdep) %*% Xdep) %*% t(Xdep) %*% Ydep
Res= Ydep- Xdep %*% bOLS
Res=Res[1:T,1]
b2=solve(t(Res) %*% Res/T) %*% (1-rho %*% rho)/(n-2)
Res=as.matrix(Res,ncol=1)
##obtain the log likelihood
LL1=lik_clo(Res,b2,n,rho)
LL1[1]
## [[1]]
## [1] -125.9855
First, save the constants obtained from evaluating the function
lik_clo
as follows:
deg=200
niter=200
AllSt=matrix(unlist(LL1[3]), ncol=1)
allctil=matrix(unlist(LL1[4]),nrow=T, ncol=(deg+1))
donde=(niter>deg)*niter+(deg>=niter)*deg
alogfac=matrix(unlist(LL1[5]),nrow=(deg+1),ncol=(donde+1))
alogfac2=matrix(unlist(LL1[6]), ncol=1)
alfac=matrix(unlist(LL1[7]), ncol=1)
repli is the number of replications. Then by averaging draws from the exact posterior distribution of the inverse volatilities, the smoothed estimates of the volatility can be obtained.
milaK=0
repli=5
keep0=matrix(0,nrow=repli, ncol=1)
for (jj in 1:repli)
{
laK=DrawK0(AllSt,allctil,alogfac, alogfac2, alfac, n, rho, b2,nproc2=2)
milaK=milaK+1/laK*(1/repli)
keep0[jj]=mean(1/laK)/b2
}
ccc=1/b2
fefo=as.vector(milaK)*ccc
## Warning in as.vector(milaK) * ccc: Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.
##obtain moving average of squared residuals
mRes=matrix(0,nrow=T,ncol=1)
Res2=Res*Res
bandi=5
for (iter in 1:T)
{ low=(iter-bandi)*(iter>bandi)+1*(iter<=bandi)
up=(iter+bandi)*(iter<=(T-bandi))+T*(iter>(T-bandi))
mRes[iter]=mean(Res2[low:up])
}
##plot the results
plot(fefo,type="l", col = "red", xlab="Time",ylab="Volatility Means")
lines(mRes, type="l", col = "blue")
legend("topright", legend = c("Stochastic Volatility", "Squared Residuals"),
col = c("red", "blue"), lty = 1, cex = 0.8)
##usage of ourgeo to evaluate a 2F1 hypergeometric function
ourgeo(1.5,1.9,1.2,0.7)
## [1] 15.93883
National Graduate Institute for Policy Studies, phd20303@grips.ac.jp↩︎
National Graduate Institute for Policy Studies, rlg@grips.ac.jp↩︎