---
title: "A Tutorial for invgamstochvol package"
author:
- Blessings Majoni^[National Graduate Institute for Policy Studies, phd20303@grips.ac.jp]
- Roberto Leon Gonzalez^[National Graduate Institute for Policy Studies, rlg@grips.ac.jp]
output: 
  html_document:
    number_sections: yes
  pdf_document: default
thanks: "This research was supported in part with\nGrant-in-Aid for Scientific Research(Kakenhi Project)
  \nfrom
  the Japan Society for the Promotion of Sciences.\n"
linkcolor: red
citecolor: green
urlcolor: blue
fontsize: 11pt
keywords: invgamstochvol
abstract: This vignette is a tutorial for a R package `invgamstochvol`. 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.
vignette: >
  %\VignetteIndexEntry{A Tutorial for invgamstochvol package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

### 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

```{r eval=TRUE}
#install.packages("invgamstochvol")
```

### Usage

```{r eval=TRUE}
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.

```{r eval=TRUE}
##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.

```{r eval=TRUE}
## 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.

```{r eval=TRUE}
## obtain the log likelihood
LL1=lik_clo(Res,b2,n,rho)
LL1[1]
```

### 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.

```{r eval = TRUE}
##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)
```

```{r eval=TRUE}
## 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)
```

```{r eval=TRUE}
##obtain the log likelihood 
LL1=lik_clo(Res,b2,n,rho)
LL1[1]
```

### Obtain smoothed estimates of volatility.

First, save the constants obtained from evaluating the function `lik_clo` as follows:

```{r eval=TRUE}
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.

```{r eval=TRUE}
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
```

```{r eval=TRUE}
##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)
```
