#_____________ MCMC routine for Gaussian AR(1) __________________
# Data statistics
xsts <- function( x )
{
T <- length(x)
xt <- x[2:T]
xt1 <- x[1:(T-1)]
xt_bar <- mean(xt)
xt1_bar <- mean(xt1)
st12 <- sum(xt1^2)
stt1 <- sum(xt*xt1)
alt1<- stt1/st12
alt2 <- ar(x,order.max=1,method="mle")
if(length(alt2$ar)==0) mle=alt1 else mle=alt2$ar
cat("MLE of slope=",format(mle,digits=3),"\n")
return(list(xt=xt,xt1=xt1,xtbar=xt_bar,xt1bar=xt1_bar,st12=st12,stt1=stt1,T=T-1,mle=mle))
}
# main routine. Gibbs sampler.
# aprior are the parameters of a Gaussian prior for the steady state
# bprior are the parameters of a Gaussian prior for the slope
# lprior are the parameters of a Gamma prior for the error precision
# x0 is the initial point
# output is a list with the chains, DIC and data
MCMC_AR1non <- function( M,dats,x0=list(a=0,b=0,l=1),aprior,bprior,lprior )
{
cat("Gibbs routine AR(1) without restrictions\n")
# initialise output
Al <- vector('numeric',M)
Be <- vector('numeric',M)
La <- vector('numeric',M)
dic <- vector('numeric',M)
# initialise chains
ca <- x0$a
cb <- x0$b
cl <- x0$l
Al[1] <- ca
Be[1] <- cb
La[1] <- cl
dats <- dats[!is.na(dats)]
# calculate statistics
S <- xsts(dats)
astar <- S$T/2+lprior[1]
dic[1] <- -(S$T*log(cl) - cl*sum((S$xt- ca - cb*S$xt1)^2))
for(i in 2:M){
pstar <- S$T*cl + aprior[2]
mstar <- (S$T*cl*(S$xtbar - cb*S$xt1bar) + prod(aprior) )/pstar
ca <- mstar + rnorm(1)/sqrt(pstar)
bstar <- lprior[2] + sum((S$xt- ca - cb *S$xt1)^2)/2
cl <- rgamma(1,astar,rate=bstar)
pstar <- cl*S$st12 + bprior[2]
mstar <- (cl*(S$stt1 - ca*S$T*S$xtbar) + prod(bprior))/pstar
cb <- mstar + rnorm(1)/sqrt(pstar)
La[i] <- cl
Be[i] <- cb
Al[i] <- ca
dic[i] <- -(S$T*log(cl) - cl*sum((S$xt- ca - cb*S$xt1)^2))
if(i%%(M/10)==0) cat(i/M*100,"% done \n")
}
DIC <- 2*mean(dic) + S$T*log(mean(La)) - mean(La)*sum((S$xt- mean(Al) - mean(Be)*S$xt1)^2)
return(list(alpha=Al,beta=Be,lambda=La,data=S,DIC=DIC))
}#_____________ MCMC routine for Gaussian AR(1) with stationarity __________________
# pdf of a re-scaled beta
rebe <- function( x,a,b )
{
p <- -(a+b-1)*log(2) - lbeta(a,b) + (a-1)*log(1 +x) + (b-1)*log(1-x)
return(p)
}
# full conditional of beta
fc_beta <- function( be,al,la,st12,stt1,st1,a,b )
{
p <- -la/2*(be^2*st12 - 2*be*(stt1-al*st1)) + (a-1)*log(be+1) + (b-1)*log(1-be)
return(p)
}
# data statistics
datasts <- function( x )
{
T <- length(x)
xt <- x[2:T]
xt1 <- x[1:(T-1)]
xt_bar <- mean(xt)
xt1_bar <- mean(xt1)
st1 <- sum(xt1)
st12 <- sum(xt1^2)
stt1 <- sum(xt*xt1)
alt1<- stt1/st12
alt2 <- ar(x,order.max=1,method="mle")
if(length(alt2$ar)==0) mle=alt1 else mle=alt2$ar
cat("MLE of slope=",format(mle,digits=3),"\n")
return(list(xt=xt,xt1=xt1,xtbar=xt_bar,xt1bar=xt1_bar,st1=st1,st12=st12,stt1=stt1,T=T-1,mle=mle))
}
# mh step for beta using a RW proposal. use v to achieve acceptance rate
move.rw_beta <- function(cb,al,la,S,a,b,v){
# if(abs(S$mle)<1){
aux <- (S$mle^2+v-1)/(2*v)
pa <- -(S$mle+1)*aux
pb <- (S$mle-1)*aux
if(pb<=0 || pa<=0){
cat("variance not allowed\n")
cat("proposal parameters are negative\n")
cat(c(pa,pb),"\n")
stop("variance too large")
}
# }
# else{
# cat(S$mle,"\n")
# pa <- 3
# pb <-3
# }
prop <- 2*rbeta(1,pa,pb)-1
mhr <- fc_beta(prop,al,la,S$st12,S$stt1,S$st1,a,b) - fc_beta(cb,al,la,S$st12,S$stt1,S$st1,a,b)
u <- log(runif(1))
alpha <- min(0,mhr)
dec <- (alpha>=u)*1
new <- prop*dec + cb*(1-dec)
return(list(value=new,acc=dec))
}
move.beta <- function(cb,al,la,S,a,b,p){
b.c <- p*(1-cb)/(1+cb)
prop <- 2*rbeta(1,p,b.c)-1
c.b <- p*(1-prop)/(1+prop)
mhr <- fc_beta(prop,al,la,S$st12,S$stt1,S$st1,a,b) - fc_beta(cb,al,la,S$st12,S$stt1,S$st1,a,b) + rebe(cb,p,c.b) - rebe(prop,p,b.c)
# print(mhr)
# print(prop)
u <- log(runif(1))
alpha <- min(0,mhr)
dec <- (alpha>=u)*1
new <- prop*dec + cb*(1-dec)
return(list(value=new,acc=dec))
}
# main routine.
# aprior are the parameters of a Gaussian prior for the steady state
# bprior are the parameters of a re-scaled Beta prior for the slope, default are (3,3)
# lprior are the parameters of a Gamma prior for the error precision
# x0 is the initial point
# v is the variance of the proposal distribution, use this to control the acceptance rate
# output is a list with the chains, DIC, acceptance rates and data
MCMC.RW_AR1stat <- function(M,dats,x0=list(a=0,b=0,l=1),aprior,bprior=c(2,2),lprior,v=0.1 )
{
cat("MCMC routine with Random Walk proposal\n")
# initialise output
Al <- vector('numeric',M)
Be <- vector('numeric',M)
La <- vector('numeric',M)
dic <- vector('numeric',M)
# initialise chains
ca <- x0$a
cb <- x0$b
if(abs(cb)>1) stop(“starting point for beta is not stationary”)
cl <- x0$l
Al[1] <- ca
Be[1] <- cb
La[1] <- cl
# calculate statistics
dats <- dats[!is.na(dats)]
S <- datasts(dats)
astar <- S$T/2+lprior[1]
dic[1] <- -(S$T*log(cl) - cl*sum((S$xt- ca - cb*S$xt1)^2))
acb <- 0
pm <- prod(aprior)
for(i in 2:M){
pstar <- S$T*cl + aprior[2]
mstar <- (S$T*cl*(S$xtbar - cb*S$xt1bar) + pm )/pstar
ca <- mstar + rnorm(1)/sqrt(pstar)
bstar <- lprior[2] + sum((S$xt- ca - cb *S$xt1)^2)/2
# if(astar)
cl <- rgamma(1,astar,rate=bstar)
aux <- move.rw_beta(cb,ca,cl,S,bprior[1],bprior[2],v)
cb <- aux$value
acb <- acb + aux$acc
La[i] <- cl
Be[i] <- cb
Al[i] <- ca
dic[i] <- -(S$T*log(cl) - cl*sum((S$xt- ca - cb*S$xt1)^2))
if(i%%(M/10)==0) cat(i/M*100,"% done \n",sep="")
}
DIC <- 2*mean(dic) + S$T*log(mean(La)) - mean(La)*sum((S$xt- mean(Al) - mean(Be)*S$xt1)^2)
cat("acc rate=",format(acb/M,digits=3),"\n")
return(list(alpha=Al,beta=Be,lambda=La,acc=acb,data=S,DIC=DIC))
}
MCMC.AR1stat <- function(M,dats,x0=list(a=0,b=0,l=1),aprior,bprior=c(2,2),lprior,p=2)
{
cat("MCMC routine with Metropolis\n")
if(p<=1) stop("Proposal parameter must be larger than 1")
# initialise output
Al <- vector('numeric',M)
Be <- vector('numeric',M)
La <- vector('numeric',M)
dic <- vector('numeric',M)
# initialise chains
ca <- x0$a
cb <- x0$b
if(abs(cb)>1) stop(“starting point for beta is not stationary”)
cl <- x0$l
Al[1] <- ca
Be[1] <- cb
La[1] <- cl
# calculate statistics
dats <- dats[!is.na(dats)]
S <- datasts(dats)
astar <- S$T/2+lprior[1]
dic[1] <- -(S$T*log(cl) - cl*sum((S$xt- ca - cb*S$xt1)^2))
acb <- 0
pm <- prod(aprior)
for(i in 2:M){
pstar <- S$T*cl + aprior[2]
mstar <- (S$T*cl*(S$xtbar - cb*S$xt1bar) + pm )/pstar
ca <- mstar + rnorm(1)/sqrt(pstar)
bstar <- lprior[2] + sum((S$xt- ca - cb *S$xt1)^2)/2
# if(astar)
cl <- rgamma(1,astar,rate=bstar)
aux <- move.beta(cb,ca,cl,S,bprior[1],bprior[2],p)
cb <- aux$value
acb <- acb + aux$acc
La[i] <- cl
Be[i] <- cb
Al[i] <- ca
dic[i] <- -(S$T*log(cl) - cl*sum((S$xt- ca - cb*S$xt1)^2))
if(i%%(M/10)==0) cat(i/M*100,"% done \n",sep="")
}
DIC <- 2*mean(dic) + S$T*log(mean(La)) - mean(La)*sum((S$xt- mean(Al) - mean(Be)*S$xt1)^2)
cat("acc rate=",format(acb/M,digits=3),"\n")
return(list(alpha=Al,beta=Be,lambda=La,acc=acb,data=S,DIC=DIC))
}# Example
load("WB_PCGR.RData")
load("CountryCodes.Rdata")
source("AR1_nonstat.R")
source("AR1_stationary.R")
# quick look at the data file
str(GDP_growthpc)
# choose a country
Cou <- "SGP"
Idx <-Ccodes$Code==Cou
Tag <- Ccodes$Region[Idx]
# Look at the series
quartz("Growth rate",width=10,height=7)
plot(GDP_growthpc[,Cou],ylab="Growth rate", xlab="Year", main=Tag,col="navyblue",lwd=2)
# _______________ Stationary model _____________________
# prior parameters and starting point
a.prior <- c(0,0.1)
l.prior <- c(2,0.6)
b.prior <- c(2,2)
ini.point <- list(a=0,b=0,l=3)
# run length, burn in and thinning
MM <- 1.05e5; bi<- 5e3; th <- 5
keep <- seq(bi+1,MM,by=th)
# main routine
MCout_st <- MCMC.RW_AR1stat(M=MM,dats=GDP_growthpc[,Cou],aprior=a.prior,bprior=b.prior,lprior=l.prior,x0=ini.point,v=0.13)
quartz("Traces Stationary",width=10,height=7)
par(mfrow=c(2,2),pty="m",mai=c(0.9,0.5,0.2,0.5))
ts.plot(MCout_st$beta,main="slope",xlab="Draw")
ts.plot(MCout_st$alpha,main="Intercept",xlab="Draw")
ts.plot(MCout_st$lambda,main="Precision",xlab="Draw")
ts.plot(MCout_st$alpha[keep]/(1-MCout_st$beta[keep]), main="Steady state",xlab="Draw")
# plot posterior distributions
quartz(paste("Priors and posteriors",Tag),width=10,height=7)
par(mfrow=c(2,2),pty="m",mai=c(0.9,0.5,0.2,0.5))
hist(MCout_st$beta[keep],nclass=80,freq=F,main="Slope",xlab=expression(beta),col="honeydew2",border="steelblue")
lines(density(MCout_st$beta[keep]),lwd=2,col="blue4")
curve(exp(rebe(x,b.prior[1],b.prior[2])),from=-1,to=1,lwd=2,col="red",lty=3,add=T)
hist(MCout_st$alpha[keep],nclass=80,freq=F,main="Intercept",xlab=expression(alpha),col="honeydew2",border="steelblue")
curve(dnorm(x,a.prior[1],1/sqrt(a.prior[2])),from=-3,to=6,lwd=2,col="red",lty=3,add=T)
lines(density(MCout_st$alpha[keep]),lwd=2,col="blue4")
hist(MCout_st$lambda[keep],nclass=80,freq=F,main="Scale",xlab=expression(lambda),col="honeydew2",border="steelblue")
curve(dgamma(x,l.prior[1],rate=l.prior[2]),from=0,to=10,lwd=2,col="red",lty=3,add=T)
lines(density(MCout_st$lambda[keep]),lwd=2,col="blue4")
hist(MCout_st$alpha[keep]/(1-MCout_st$beta[keep]),nclass=90,freq=F,main="Steady state", xlab=expression(rho), col="honeydew2", border="steelblue",xlim=c(0,10),axes=F)#, ylim=c(0,13),axes=F,xlim=c(-0.15,0.5))
axis(1, at=seq(-6,12,by=1))
axis(2, at=seq(0,2,by=0.25))
lines(density(MCout_st$alpha[keep]/(1-MCout_st$beta[keep])),lwd=2,col="blue4")
# _______________ Non-Stationary model _____________________
# prior parameters and starting point
a.prior <- c(0,0.1)
b.prior <- c(0,0.1)
l.prior <- c(2,0.6)
ini.point <- list(a=0,b=5,l=1)
# run length, burn in and thinning
MM <- 1.05e5; bi<- 5e3; th <- 5
keep <- seq(bi,MM,by=th)
# main routine
MCout_ns <- MCMC_AR1non(M=MM,dats=GDP_growthpc[,Cou],aprior=a.prior,bprior=b.prior,lprior=l.prior,x0=ini.point)
quartz(paste("Traces NS",Tag),width=10,height=7)
par(mfrow=c(2,2),pty="m",mai=c(0.9,0.5,0.2,0.5))
ts.plot(MCout_ns$beta,main="slope",xlab="Draw")
ts.plot(MCout_ns$alpha,main="Intercept",xlab="Draw")
ts.plot(MCout_ns$lambda,main="Precision",xlab="Draw")
# plot ppsterior distributions
quartz(paste("Priors and posteriors NS",Tag),width=10,height=7)
par(mfrow=c(2,2),pty="m",mai=c(0.9,0.5,0.2,0.5))
hist(MCout_ns$beta[keep],nclass=80,freq=F,main="Slope",xlab=expression(beta),col="lightyellow2",border="steelblue")
lines(density(MCout_ns$beta[keep]),lwd=2,col="blue4")
curve(dnorm(x,b.prior[1],1/sqrt(b.prior[2])),from=-4,to=8,lwd=2,col="red",lty=3,add=T)
hist(MCout_ns$alpha[keep],nclass=80,freq=F,main="Intercept",xlab=expression(alpha),col="lightyellow2",border="steelblue")
lines(density(MCout_ns$alpha[keep]),lwd=2,col="blue4")
curve(dnorm(x,a.prior[1],1/sqrt(a.prior[2])),from=-4,to=8,lwd=2,col="red",lty=3,add=T)
hist(MCout_ns$lambda[keep],nclass=80,freq=F,main="Scale",xlab=expression(lambda),col="lightyellow2",border="steelblue")
curve(dgamma(x,l.prior[1],rate=l.prior[2]),from=0,to=2,lwd=2,col="red",lty=3,add=T)
lines(density(MCout_ns$lambda[keep]),lwd=2,col="blue4")