##################################################################################### # # Application 1: Dynamic Beta Regression # # Lopes and Tsay (2011) # Particle Filters and Bayesian Inference in Financial Econometrics. # Journal of Forecasting. # ##################################################################################### # # Data: Brazilian monthly unemployment rate from March 2002 to December 2009. # Source: Brazilian Institute for Geography and Statistics # UFL: http://www.sidra.ibge.gov.br/bda/pesquisas/pme/default.asp#dead # ###################################################################################### # # Author : Hedibert Freitas Lopes # Associate Professor of Econometrics and Statistics # The University of Chicago Booth School of Business # 5807 South Woodlawn Avenue # Chicago, Illinois, 60637 # Email : hlopes@ChicagoBooth.edu # URL: http://faculty.chicagobooth.edu/hedibert.lopes # ##################################################################################### rm(list=ls()) PL = function(y,m0,C0,a0,b0,c0,d0,N){ n = length(y) x = rnorm(N,m0,sqrt(C0)) sig2 = 1/rgamma(N,a0,b0) tau2 = 1/rgamma(N,c0,d0) S = matrix(c(a0,b0,c0,d0),N,4,byrow=T) xs = matrix(0,n,N) sig2s = matrix(0,n,N) tau2s = matrix(0,n,N) like = rep(0,n) for (t in 1:n){ A = tau2/(tau2+sig2) w = dnorm(y[t],x,sqrt(sig2+tau2)) like[t] = sum(w) k = sample(1:N,size=N,prob=w,replace=T) x1 = x[k] A = A[k] sig2 = sig2[k] S = S[k,] x = rnorm(N,A*y[t]+(1-A)*x1,sqrt(A*sig2)) S[,1] = S[,1] + 0.5 S[,2] = S[,2] + 0.5*(y[t]-x)^2 S[,3] = S[,3] + 0.5 S[,4] = S[,4] + 0.5*(x-x1)^2 sig2 = 1/rgamma(N,S[,1],S[,2]) tau2 = 1/rgamma(N,S[,3],S[,4]) xs[t,] = x sig2s[t,] = sig2 tau2s[t,] = tau2 } return(list(xs=xs,sigs=sqrt(sig2s),taus=sqrt(tau2s),like=like)) } LW = function(y,m0,C0,a0,b0,c0,d0,a,N){ n = length(y) h = sqrt(1-a^2) beta = rnorm(N,m0,sqrt(C0)) phi = 1/rgamma(N,a0,b0) tau2 = 1/rgamma(N,c0,d0) pars = cbind(log(phi),log(tau2)) draws = array(0,c(n,N,3)) like = rep(0,n) for (t in 1:n){ mu = 1/(1+exp(-beta)) w = dbeta(y[t],exp(pars[,1])*mu,exp(pars[,1])*(1-mu)) like[t] = sum(w) mean = apply(pars,2,mean) vari = var(pars) m1 = a*pars+(1-a)*mean w = dbeta(y[t],exp(m1[,1])*mu,exp(m1[,1])*(1-mu),log=TRUE) w = exp(w-max(w)) k = sample(1:N,replace=TRUE,size=N,prob=w) pars1 = pars[k,]+matrix(rnorm(N*2),N,2)%*%chol(vari)*h beta1 = rnorm(N,beta[k],exp(pars1[,2]/2)) mu1 = 1/(1+exp(-beta1)) w1 = dbeta(y[t],exp(pars1[,1])*mu1,exp(pars1[,1])*(1-mu1),log=TRUE) - log(w[k]) w1 = exp(w1-max(w1)) k = sample(1:N,replace=TRUE,size=N,prob=w1) pars = pars1[k,] beta = beta1[k] draws[t,,] = cbind(sqrt(1/(1+exp(pars[,1]))),1/(1+exp(-beta)),exp(pars[,2]/2)) } return(list(draws=draws,like=like)) } y = c(2.566,2.488,2.353,2.298,2.377,2.381,2.337,2.283,2.238,2.130,2.306,2.399,2.527,2.606,2.716,2.750, 2.697,2.783,2.795,2.780,2.622,2.314,2.455,2.535,2.737,2.823,2.634,2.524,2.421,2.471,2.379,2.282,2.330, 2.080,2.184,2.291,2.352,2.341,2.221,2.029,2.041,2.045,2.118,2.110,2.109,1.823,2.016,2.201,2.281,2.266, 2.229,2.306,2.394,2.384,2.257,2.209,2.150,1.864,2.063,2.195,2.281,2.272,2.274,2.187,2.137,2.178,2.058, 1.986,1.886,1.681,1.806,1.970,1.952,1.948,1.802,1.807,1.867,1.752,1.777,1.743,1.760,1.567,1.890,1.941, 2.082,2.046,2.036,1.867,1.854,1.889,1.799,1.753,1.714,1.592)/20 n = length(y) month = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec") month = c(month[3:12],rep(month,7)) year = matrix(matrix(2002:2009,12,8,byrow=TRUE),12*8,1)[3:(n+2),1] data = cbind(month,year,y) ind = c(1,seq(11,n,by=12),n) lab = c("mar/03","jan/03","jan/04","jan/05","jan/06","jan/07","jan/08","jan/09","dec/09") plot(y,axes=F,xlab="Months",ylab="Unemployment",pch=16,type="b") axis(2,at=seq(0.06,0.14,by=0.02)) axis(1,at=ind,lab=lab) box() set.seed(4654657) N = 10000 m0 = 0.1 C0 = 100 a0 = 2.1 b0 = (a0+1)*0.00001 c0 = 2.1 d0 = (a0+1)*0.00001 pl = PL(y,m0,C0,a0,b0,c0,d0,N) qs = array(0,c(3,n,3)) qs[1,,]= t(apply(pl$xs,1,quantile,c(0.025,0.5,0.975))) qs[2,,]= t(apply(pl$sigs,1,quantile,c(0.025,0.5,0.975))) qs[3,,]= t(apply(pl$taus,1,quantile,c(0.025,0.5,0.975))) m0 = log(y[1]/(1-y[1])) C0 = 0.1 a0 = 2.1 b0 = (a0+1)*15000 c0 = 2.1 d0 = (c0+1)*0.05 a = 0.95 lw = LW(y,m0,C0,a0,b0,c0,d0,a,N) lw1 = lw$draws qs1 = array(0,c(3,n,3)) qs1[1,,] = t(apply(lw1[,,1],1,quantile,c(0.025,0.5,0.975))) qs1[2,,] = t(apply(lw1[,,2],1,quantile,c(0.025,0.5,0.975))) qs1[3,,] = t(apply(lw1[,,3],1,quantile,c(0.025,0.5,0.975))) pdf(file="models.pdf",width=8,height=10) par(mfrow=c(3,2)) plot(qs[1,,1],axes=F,xlab="Months",ylab="",main="(a)",type="l",ylim=range(qs[1,,],qs1[2,,]),col=1,lwd=2) axis(2,at=seq(0.06,0.14,by=0.02));axis(1,at=ind,lab=lab);box() for (i in 2:3) lines(qs[1,,i],col=1,lwd=2) for (i in 1:3) lines(qs1[2,,i],col=grey(0.5),lwd=2) legend(40,0.15,legend=c("Local level","Dynamic Beta"),cex=1.5,col=c(1,grey(0.5)),lwd=c(2,2),lty=c(1,1),bty="n") BF = cumsum(lw$like)/cumsum(pl$like) plot(BF,axes=F,xlab="Months",ylab="Bayes factor",main="(b)",type="l",pch=16,ylim=range(BF[2:n],1)) axis(2);axis(1,at=ind,lab=lab);box() abline(h=1,lty=2) plot(qs[2,,2],axes=F,xlab="Months",main="(c)",ylab=expression(sigma),pch=16,ylim=range(qs[2,,]),type="l") axis(2);axis(1,at=ind,lab=lab);box() for (i in c(1,3)) lines(qs[2,,i],col=grey(0.5)) plot(qs[3,,2],axes=F,xlab="Months",main="(d)",ylab=expression(tau),pch=16,ylim=range(qs[3,,]),type="l") axis(2);axis(1,at=ind,lab=lab);box() for (i in c(1,3)) lines(qs[3,,i],col=grey(0.5)) plot(qs1[1,,2],axes=F,xlab="Months",main="(e)",ylab=expression((1/(1+phi))^(1/2)),pch=16,ylim=range(qs1[1,,]),type="l") axis(2);axis(1,at=ind,lab=lab);box() for (i in c(1,3)) lines(qs1[1,,i],col=grey(0.5)) plot(qs1[3,,2],axes=F,xlab="Months",main="(f)",ylab=expression(W^(1/2)),pch=16,ylim=range(qs1[3,,]),type="l") axis(2);axis(1,at=ind,lab=lab);box() for (i in c(1,3)) lines(qs1[3,,i],col=grey(0.5)) dev.off() # MONTE CARLO STUDY # ----------------- N = 10000 rep = 20 qs = array(0,c(rep,2,2,n,3)) BF = matrix(0,rep,n) set.seed(13164) for (i in 1:rep){ m0 = 0.1 C0 = 100 a0 = 2.1 b0 = (a0+1)*0.00001 c0 = 2.1 d0 = (a0+1)*0.00001 pl = PL(y,m0,C0,a0,b0,c0,d0,N) qs[i,1,1,,]= t(apply(pl$sigs,1,quantile,c(0.025,0.5,0.975))) qs[i,1,2,,]= t(apply(pl$taus,1,quantile,c(0.025,0.5,0.975))) m0 = log(y[1]/(1-y[1])) C0 = 0.1 a0 = 2.1 b0 = (a0+1)*15000 c0 = 2.1 d0 = (c0+1)*0.05 a = 0.95 lw = LW(y,m0,C0,a0,b0,c0,d0,a,N) lw1 = lw$draws qs[i,2,1,,] = t(apply(lw1[,,1],1,quantile,c(0.025,0.5,0.975))) qs[i,2,2,,] = t(apply(lw1[,,3],1,quantile,c(0.025,0.5,0.975))) BF[i,] = cumsum(lw$like)/cumsum(pl$like) } pdf(file="MCstudy.pdf",width=8,height=10) cols = c(grey(0.5),1,grey(0.5)) par(mfrow=c(2,2)) plot(qs[1,1,1,,1],xlab="",ylab="",main=expression(sigma),ylim=range(qs[,1,1,,]),type="l",axe=F,col=cols[1]) axis(2);axis(1,at=ind,lab=lab);box() for (i in 2:3) lines(qs[1,1,1,,i],col=cols[i]) for (j in 2:rep) for (i in 1:3) lines(qs[j,1,1,,i],col=cols[i]) plot(qs[1,1,2,,1],xlab="",ylab="",main=expression(tau),ylim=range(qs[,1,2,,]),type="l",axe=F,col=cols[1]) axis(2);axis(1,at=ind,lab=lab);box() for (i in 2:3) lines(qs[1,1,2,,i],col=cols[i]) for (j in 2:rep) for (i in 1:3) lines(qs[j,1,2,,i],col=cols[i]) plot(qs[1,2,1,,1],xlab="",ylab="",main=expression((1/(1+phi))^(1/2)),ylim=range(qs[,2,1,,]),type="l",axe=F,col=cols[1]) axis(2);axis(1,at=ind,lab=lab);box() for (i in 2:3) lines(qs[1,2,1,,i],col=cols[i]) for (j in 2:rep) for (i in 1:3) lines(qs[j,2,1,,i],col=cols[i]) plot(qs[1,2,2,,1],xlab="",ylab="",main=expression(W^(1/2)),ylim=range(qs[,2,2,,]),type="l",axe=F,col=cols[1]) axis(2);axis(1,at=ind,lab=lab);box() for (i in 2:3) lines(qs[1,2,2,,i],col=cols[i]) for (j in 2:rep) for (i in 1:3) lines(qs[j,2,2,,i],col=cols[i]) dev.off() pdf(file="BF.pdf",width=5,height=4) par(mfrow=c(1,1)) plot(BF[1,2:n],xlab="",ylab="Bayes factor",main="",ylim=range(BF[,2:n]),type="l",axe=F,col=grey(0.5)) axis(2);axis(1,at=ind,lab=lab);box() for (j in 2:rep) lines(BF[j,2:n],col=grey(0.5)) dev.off()