############################################################################################### # # Example ii. LOG STOCHASTIC VOLATILITY AR(1) MODEL # # AUXILIARY PARTICLE FILTER + PARAMETER LEARNING # ############################################################################################### # # Author : Hedibert Freitas Lopes # Graduate School of Business # University of Chicago # 5807 South Woodlawn Avenue # Chicago, Illinois, 60637 # Email : hlopes@ChicagoGSB.edu # ############################################################################################### pri = function(theta,theta0,sV0){dnorm(theta,theta0,sV0)} likelihood = function(y,theta){dnorm(y,0,exp(theta/2))} rlike = function(theta){rnorm(1,0,exp(theta/2))} post = function(y,theta,theta0,V0){pri(theta,theta0,sV0)*likelihood(y,theta)} quant025 = function(x){quantile(x,0.025)} quant975 = function(x){quantile(x,0.975)} # Simulating the data # ------------------- set.seed(12345) n = 500 alpha = -0.0031 beta = 0.9951 tau2 = 0.0074 alpha.true = alpha beta.true = beta tau2.true = tau2 tau = sqrt(tau2) y = rep(0,n) theta = rep(0,n) theta[1] = alpha/(1-beta) y[1] = rlike(theta[1]) for (t in 2:n){ theta[t] = rnorm(1,alpha+beta*theta[t-1],tau) y[t] = rlike(theta[t]) } # Data and prior hyperparameters # ------------------------------ theta0 = 0.0 V0 = 0.1 sV0 = sqrt(V0) ealpha = alpha valpha = 0.01 ephi = beta vphi = 0.01 nu = 3 lambda = tau2 # Auxiliary particle filter # ------------------------- set.seed(8642) N = 1000 thetass = NULL parss = array(0,c(N,3,n)) ws = NULL thetas = rnorm(N,theta0,sV0) pars = cbind(rnorm(N,ealpha,sqrt(valpha)),rnorm(N,ephi,sqrt(vphi)),log(1/rgamma(N,nu/2,nu*lambda/2))) delta = 0.75 h2 = 1-((3*delta-1)/(2*delta))^2 a = sqrt(1-h2) ESS = NULL par(mfrow=c(1,1)) for (t in 1:n){ mpar = apply(pars,2,mean) vpar = var(pars) ms = a*pars+(1-a)*matrix(mpar,N,3,byrow=T) mus = pars[,1]+pars[,2]*thetas weight = likelihood(y[t],mus) k = sample(1:N,size=N,replace=T,prob=weight) ms1 = ms[k,] + matrix(rnorm(3*N),N,3)%*%chol(h2*vpar) thetat = rnorm(N,ms1[,1]+ms1[,2]*thetas[k],exp(ms1[,3]/2)) w = likelihood(y[t],thetat)/likelihood(y[t],mus[k]) w = w/sum(w) ind = sample(1:N,size=N,replace=T,prob=w) thetas = thetat[ind] pars = ms1[ind,] thetass = rbind(thetass,thetas) parss[,,t] = pars ws = rbind(ws,w) cv2 = var(w)/(mean(w)^2) ESS = c(ESS,N/(1+cv2)) ts.plot(ESS,xlim=c(1,n)) } mvar.smc = apply(exp(thetass),1,mean) lvar.smc = apply(exp(thetass),1,quant025) uvar.smc = apply(exp(thetass),1,quant975) malpha = apply(parss[,1,],2,mean) lalpha = apply(parss[,1,],2,quant025) ualpha = apply(parss[,1,],2,quant975) mbeta = apply(parss[,2,],2,mean) lbeta = apply(parss[,2,],2,quant025) ubeta = apply(parss[,2,],2,quant975) mtau2 = apply(exp(parss[,3,]),2,mean) ltau2 = apply(exp(parss[,3,]),2,quant025) utau2 = apply(exp(parss[,3,]),2,quant975) # Graphical analysis # ------------------ par(mfrow=c(2,2)) ts.plot(mvar.smc,ylim=c(0,max(uvar.smc)),xlab="time",ylab="",main="") title(paste("SMC: N=",N,sep="")) lines(mvar.smc,col=2,lwd=1.25) lines(lvar.smc,col=2,lwd=1.25) lines(uvar.smc,col=2,lwd=1.25) lines(0.2*abs(y),type="h",col=1,lwd=2) legend(200,12,legend=c("Seq.MCMC","APF"),lwd=c(1.25,1.25),col=c(4,2),bty="n") ts.plot(malpha,ylim=range(lalpha,ualpha),ylab="",main=expression(alpha)) lines(lalpha,col=2,lwd=1.25) lines(malpha,col=2,lwd=1.25) lines(ualpha,col=2,lwd=1.25) abline(h=alpha.true,col=1,lwd=1.25) ts.plot(mbeta,ylim=range(lbeta,ubeta),ylab="",main=expression(beta)) lines(lbeta,col=2,lwd=1.25) lines(mbeta,col=2,lwd=1.25) lines(ubeta,col=2,lwd=1.25) abline(h=beta.true,col=1,lwd=1.25) ts.plot(mtau2,ylim=range(ltau2,utau2),ylab="",main=expression(tau^2)) lines(ltau2,col=2,lwd=1.25) lines(mtau2,col=2,lwd=1.25) lines(utau2,col=2,lwd=1.25) abline(h=tau2.true,col=1,lwd=1.25)