############################################################################################ # # AR(1) PLUS NOISE MODEL - SUBROUTINES # # Observation equation: y(t)|x(t),theta ~ N(x(t),sig2) # Evolution equation: x(t)|x(t-1),theta ~ N(alpha+beta*x(t-1),tau2) # # where theta=(alpha,beta), x(0) ~ N(m0,C0) and theta ~ N(a0,tau2*A0), for known a0 and # A0=diag(A0[1,1],A0[2,2]). The pair (sig2,tau2) is knwon throughout this exercise. # ############################################################################################ # # DLM.SIM - simulate data # # DLM - Means and variances of p(x_t|y^t,theta) and p(x_t|y^n,theta) plus p(y^n|theta). # # FFBS - Sampling from p(x^n|y^n,theta). # # OAPF - Sampling from p(x_t|y^t,theta) via fully adapted APF. # # PS - Sampling from p(x_t|y^n,theta) via backward particle sampling. # # MCMC - Sampling from p(x^n,theta|y^n) via Gibbs sampler: # p(x^n|theta,y^n) via FFBS. # p(theta|x^n,y^n) via Bayesian linear regression updates. # # MCMC1 - Sampling from p(x^n,theta|y^n): # p(theta|y^n) via random-walk Metropolis (analytically integrating out x^n), # p(x^n|y^n,theta) via FFBS. # # MCMC2 - Sampling from p(x^n,theta|y^n): # p(theta|y^n) via random-walk Metropolis (integrating out x^n via OAPF), # p(x^n|y^n,theta) via FFBS. # # LW - Sampling from p(x_t,theta|y^t) via LW filter. # # PL - Sampling from p(x_t,theta|y^t) via particle learning. # # STORVIK - Sampling from p(x_t,theta|y^t) via Storvik's filter. # # PS1 - Sampling from p(x_t|y^n) via backward particle sampling. # ############################################################################################ # # 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/research/ # ############################################################################################ DLM.SIM = function(n,tau,sig,x0){ w = rnorm(n,0,tau) v = rnorm(n,0,sig) x = rep(0,n) y = rep(0,n) x[1] = alpha+beta*x0 + w[1] y[1] = x[1] + v[1] for (t in 2:n){ x[t] = alpha+beta*x[t-1] + w[t] y[t] = x[t] + v[t] } return(list(x=x,y=y)) } DLM = function(y,alpha,beta,sig2,tau2,m0,C0){ n = length(y) mf = rep(0,n) Cf = rep(0,n) mb = rep(0,n) Cb = rep(0,n) a = rep(0,n) R = rep(0,n) B = rep(0,n-1) R[1] = beta^2*C0 + tau2 a[1] = alpha + beta*m0 Q = R[1] + sig2 l = dnorm(y[1],a[1],sqrt(Q),log=TRUE) A = R[1]/Q mf[1] = (1-A)*a[1]+A*y[1] Cf[1] = A*sig2 for (t in 2:n){ R[t] = beta^2*Cf[t-1] + tau2 B[t-1] = Cf[t-1]*beta/R[t] a[t] = alpha + beta*mf[t-1] Q = R[t] + sig2 l = l + dnorm(y[t],a[t],sqrt(Q),log=TRUE) A = R[t]/Q mf[t] = (1-A)*a[t]+A*y[t] Cf[t] = A*sig2 } mb[n] = mf[n] Cb[n] = Cf[n] for (t in (n-1):1){ mb[t] = mf[t] + B[t]*(mb[t+1]-a[t+1]) Cb[t] = Cf[t] + B[t]^2*(Cb[t+1]-R[t+1]) } return(list(mf=mf,Cf=Cf,mb=mb,Cb=Cb,l=l)) } FFBS = function(y,alpha,beta,sig2,tau2,m0,C0,M){ n = length(y) mf = rep(0,n) Cf = rep(0,n) a = rep(0,n) R = rep(0,n) B = rep(0,n-1) R[1] = beta^2*C0 + tau2 a[1] = alpha + beta*m0 Q = R[1] + sig2 l = dnorm(y[1],a[1],sqrt(Q),log=TRUE) A = R[1]/Q mf[1] = (1-A)*a[1]+A*y[1] Cf[1] = A*sig2 for (t in 2:n){ R[t] = beta^2*Cf[t-1] + tau2 B[t-1] = Cf[t-1]*beta/R[t] a[t] = alpha + beta*mf[t-1] Q = R[t] + sig2 l = l + dnorm(y[t],a[t],sqrt(Q),log=TRUE) A = R[t]/Q mf[t] = (1-A)*a[t]+A*y[t] Cf[t] = A*sig2 } x = matrix(0,M,n) x[,n] = rnorm(M,mf[n],sqrt(Cf[n])) for (t in (n-1):1){ mb = mf[t] + B[t]*(x[,t+1]-a[t+1]) Cb = Cf[t] - B[t]^2*R[t+1] x[,t] = rnorm(M,mb,sqrt(Cb)) } return(x) } OAPF = function(y,alpha,beta,sig2,tau2,m0,C0,M){ n = length(y) w2 = sig2+tau2 A = tau2/w2 s = sqrt(w2) d = sqrt(A*sig2) l = rep(0,n) x = rnorm(M,m0,sqrt(C0)) xs = matrix(0,n,M) for (t in 1:n){ a = alpha+beta*x w = dnorm(y[t],a,s) k = sample(1:M,size=M,prob=w,replace=T) m = A*y[t]+(1-A)*a[k] x = rnorm(M,m,d) l[t] = mean(w) xs[t,] = x } return(list(xs=xs,l=l)) } PS = function(alpha,beta,tau2,xs){ N = ncol(xs) xss = matrix(0,n,N) tau = sqrt(tau2) for (i in 1:N){ xss[n,i] = xs[n,i] for (t in (n-1):1){ w = dnorm(xss[t+1,i],alpha+beta*xs[t,],tau) k = sample(1:N,size=1,prob=w) xss[t,i] = xs[t,k] } } return(xss) } MCMC = function(y,alpha,beta,sig2,tau2,a0,A0,m0,C0,burn,niter){ n = length(y) iA = solve(A0) iAa = iA%*%a0 draw = c(alpha,beta) draws = matrix(0,niter,n+2) for (iter in 1:niter){ #print(iter) x = FFBS(y,draw[1],draw[2],sig2,tau2,m0,C0,1)[1,] V = 1/(1/C0+beta^2/tau2) m = V*(m0/C0+beta*(x[1]-alpha)/tau2) x0 = rnorm(1,m,sqrt(V)) Y = x X = cbind(1,c(x0,x[1:(n-1)])) V = solve(iA + t(X)%*%X) m = V%*%(iAa + t(X)%*%Y) V = t(chol(V)) draw = m + tau*V%*%rnorm(2) draws[iter,] = c(x,draw) } return(draws[(burn+1):niter,]) } MCMC1 = function(y,alpha,beta,sig2,tau2,a0,A0,m0,C0,burn,niter,sd){ alphaold = alpha betaold = beta den = DLM(y,alphaold,betaold,sig2,tau2,m0,C0)$l+ dnorm(alphaold,a0[1],sqrt(tau2*A0[1,1]),log=TRUE)+ dnorm(betaold,a0[2],sqrt(tau2*A0[2,2]),log=TRUE) draws1 = matrix(0,niter,n+2) for (iter in 1:niter){ print(iter) alphanew = alphaold + rnorm(1,0,sd) betanew = betaold + rnorm(1,0,sd) num = DLM(y,alphanew,betanew,sig2,tau2,m0,C0)$l+ dnorm(alphanew,a0[1],sqrt(tau2*A0[1,1]),log=TRUE)+ dnorm(betanew,a0[2],sqrt(tau2*A0[2,2]),log=TRUE) logaccep = num-den if (log(runif(1))