A Tutorial for invgamstochvol package

Introduction

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.

Installation

#install.packages("invgamstochvol")

Usage

library(invgamstochvol)

Example using simulated data

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)

Obtain the residuals

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)

Obtain the likelihood

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] -44.15683

Example using real data

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

Obtain smoothed estimates of volatility.

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)

Obtain the smoothed estimates of the volatility

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, [email protected]↩︎
National Graduate Institute for Policy Studies, [email protected]↩︎