set.seed(111)
nobs = 80; # n is often larger than 1000 in practice.
mu = 0; phi = 0.97; sigma_eta = 0.3; rho = -0.3;
h = 0; Y = c();
for(i in 1:nobs){
eps = rnorm(1, 0, 1)
eta = rho*sigma_eta*eps + sigma_eta*sqrt(1-rho^2)*rnorm(1, 0, 1)
y = eps * exp(0.5*h)
h = mu + phi * (h-mu) + eta
Y = append(Y, y)
}
# This is a toy example. Increase nsim and nburn
# until the convergence of MCMC in practice.
nsim = 300; nburn = 100;
vhyper = c(0.0,1000,1.0,1.0,1.0,1.0,0.01,0.01)
out = asv_mcmc(Y, nsim, nburn, vhyper)
vmu = out[[1]]; vphi = out[[2]]; vsigma_eta = out[[3]]; vrho = out[[4]];mh = out[[5]];
mu = mean(vmu); phi = mean(vphi); sigma_eta = mean(vsigma_eta);
rho = mean(vrho);
#
h = mh[nsim,]
theta = c(vmu[nsim],vphi[nsim],vsigma_eta[nsim],vrho[nsim])
theta_star = c(mu, phi, sigma_eta, rho)
# Increase iM in practice (such as iI = 5000, iM =5000).
result = asv_logML(h, theta, theta_star, Y, 100, 100, vHyper = vhyper)
result1 = matrix(0, 4, 2)
result1[,1] =result[[1]]
result1[,2] =result[[2]]
colnames(result1) = c("Estimate", "Std Err")
rownames(result1) = c("Log marginal lik", "Log likelihood", "Log prior", "Log posterior")
print(result1, digit=4)
Run the code above in your browser using DataLab