Log chi-square distribution

Another SIR example

1 Sampling from the log-\(\chi^2_1\)

Suppose that our interest is in drawing from \(\theta\) following a log-\(\chi^2_1\): \[ p(\theta) = \frac{1}{\sqrt 2\Gamma(0.5)}\exp\left\{0.5\theta-0.5e^\theta\right\}, \qquad \theta \in \Re. \] This distribution appears frequently in financial econometrics applications.

target = function(theta){
  1/sqrt(2)/gamma(1/2)*exp(0.5*theta-0.5*exp(theta))
}
theta = seq(-15,5,length=1000)

par(mfrow=c(1,1))
plot(theta,target(theta),xlab=expression(theta),ylab="Density",
     main="Target",type="l")

2 SIR in action

2.1 We have to pick a proposal density

Let us compare the performance of the SIR based on two alternative proposal densities: \[ q_1(\theta) \equiv Normal(0,5^2) \ \ \ \mbox{and} \ \ \ q_2(\theta) \equiv \mbox{Uniform}(-15,15). \]

q1 = function(theta){
  dnorm(theta,0,5)
}
q2 = function(theta){
  rep(1/30,length(theta))
}

par(mfrow=c(1,1))
plot(theta,target(theta),xlab=expression(theta),ylab="Density",main="Target vs proposals",type="l")
lines(theta,q1(theta),col=2,lwd=3)
lines(theta,q2(theta),col=3,lwd=3)
legend("topleft",legend=c("Target","Gaussian proposal","Uniform proposal"),col=1:3,lty=1)

2.2 Sampling from proposals

M = 1000000
theta.q1 = rnorm(M,0,5)
theta.q2 = runif(M,-15,15)

2.3 Computing acceptance ratios

w1 = target(theta.q1)/q1(theta.q1)
w2 = target(theta.q2)/q2(theta.q2)

2.4 Resampling

theta.t1 = sample(theta.q1,size=M/2,replace=TRUE,prob=w1)
theta.t2 = sample(theta.q2,size=M/2,replace=TRUE,prob=w2)

2.5 Histograms against exact density

par(mfrow=c(2,2))
hist(theta.q1,prob=TRUE,xlab=expression(theta),main="Gaussian proposal",
     xlim=c(-15,15))
hist(theta.t1,prob=TRUE,xlab=expression(theta),main="Posterior",
     xlim=c(-15,15))
lines(theta,target(theta),col=2,lwd=3)
hist(theta.q2,prob=TRUE,xlab=expression(theta),main="Uniform proposal",
     xlim=c(-15,15))
hist(theta.t2,prob=TRUE,xlab=expression(theta),main="Posterior",
     xlim=c(-15,15))
lines(theta,target(theta),col=2,lwd=3)

2.6 Computing \(Pr(\theta<-5)\)

# Exact value (via Riemann integration)
thetas = seq(-20,-5,length=10000)
h = thetas[2]-thetas[1]
I = sum(h*target(thetas))

I1 = mean(theta.t1<(-5))
I2 = mean(theta.t2<(-5))
c(I,I1,I2)

## [1] 0.06540915 0.06490600 0.06541400

3 Repeating SIR 200 times

set.seed(31416)
M   = 100000
Rep = 100
Is = matrix(0,Rep,2)
for (rep in 1:Rep){
theta.q1 = rnorm(M,0,5)
theta.q2 = runif(M,-15,15)
w1 = target(theta.q1)/q1(theta.q1)
w2 = target(theta.q2)/q2(theta.q2)
theta.t1 = sample(theta.q1,size=M/2,replace=TRUE,prob=w1)
theta.t2 = sample(theta.q2,size=M/2,replace=TRUE,prob=w2)
I1 = mean(theta.t1<(-5))
I2 = mean(theta.t2<(-5))
Is[rep,] = c(I1,I2)
}
par(mfrow=c(1,1))
boxplot(Is,names=c("Gaussian proposal","Uniform proposal"))
abline(h=I,col=2,lwd=3)

3.1 Computing root MSE

sqrt(var(Is[,1]-I))

## [1] 0.00124749

sqrt(var(Is[,2]-I))

## [1] 0.001491914

sqrt(var(Is[,1]-I))/sqrt(var(Is[,2]-I))

## [1] 0.836168