程序代写代做 chain #_____________ MCMC routine for Gaussian AR(1) __________________

#_____________ 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")