R Code
Click [+] to expand or collapse section. Click [-] to collapse entire page. Or [+ Expand the Entire Page +]
There is an R package for the text called nltsa that is available on GitHub. To install the package, issue the following two lines:
install.packages("devtools")
devtools::install_github("nickpoison/nltsa")
[+] Chapter 1
Example 1.36
library(astsa) # load astsa package sunspot = sqrt(sunspot.year) # transform sunspot data plot(sunspot, xlab="year") # plot data acf2(sunspot) # not shown ar(sunspot) # uses AIC to find the best fit with Yule-Walker estimates arima(sunspot, order=c(9,0,0)) # final model fit via MLE
Example 1.37
sunspot = sqrt(sunspot.year) # if you didn't save it from previous example
sun.ar = spec.ar(sunspot, plot=FALSE) # parametric
sun.np = spectrum(sunspot, spans=c(5,5), plot=FALSE) # nonparametric
plot(sun.ar$freq, sun.ar$spec, type="l", log="y", ylab="spectrum", xlab="frequency")
lines(sun.np$freq, sun.np$spec, lty=2)
legend("topright", c("parametric","nonparametric"), lty=1:2)
[+] Chapter 2
Example 2.8
library(astsa) # Generate data set.seed(1) num = 10 w = rnorm(num+1) v = rnorm(num) mu = cumsum(w) # states y = mu[-1] + v # observations # Smooth mu0 = 0; sigma0 = 1; phi = 1; cQ = 1; cR = 1 (ks = Ksmooth0(num, y, 1, mu0, sigma0, phi, cQ, cR)) # Plot smoother plot.ts(mu[-1], type="p", ylim=c(-3,5), main="Smoother") lines(ks$xs) lines(ks$xs+2*sqrt(ks$Ps), lty=2, col=4) lines(ks$xs-2*sqrt(ks$Ps), lty=2, col=4)
Example 2.11
library(astsa)
# Generate Data
set.seed(999); num = 100
x = arima.sim(n=num+1, list(ar = .8, sd=1))
y = ts(x[-1] + rnorm(num,0,1))
init.par = c(phi=.1, sigw=1, sigv=1) # initial parameters
# Function to evaluate the likelihood
Linn = function(para){
phi = para[1]; sigw = para[2]; sigv = para[3]
Sigma0 = (sigw^2)/(1-phi^2); Sigma0[Sigma0<0]=0
kf = Kfilter0(num, y, 1, mu0=0, Sigma0, phi, sigw, sigv)
return(kf$like)
}
# Estimation
(est = optim(init.par,Linn,gr=NULL,method="BFGS",hessian=TRUE,control=list(trace=1,REPORT=1)))
SE = sqrt(diag(solve(est$hessian)))
cbind(estimate=c(est$par[1],est$par[2],est$par[3]), SE)
Example 2.12
library(astsa)
library(nlme)
# Generate data (same as previous example)
set.seed(999); num = 100
x = arima.sim(n=num+1, list(ar = .8, sd=1))
y = ts(x[-1] + rnorm(num,0,1))
# EM procedure - output not shown
(em = EM0(num, y, 1, mu0=0, Sigma0=2, Phi=.1, cQ=1, cR=1, 75, .00001))
# Standard Errors (this uses nlme)
mu0 = em$mu0; Sigma0 = em$Sigma0
para = c(em$Phi, sqrt(em$Q), sqrt(em$R))
# evaluate likelihood at estimates
Linn = function(para){
kf = Kfilter0(num, y, 1, mu0, Sigma0, para[1], para[2], para[3])
return(kf$like)
}
# get SEs
emhess = fdHess(para, function(para) Linn(para))
SE = sqrt(diag(solve(emhess$Hessian)))
# Display Summary of Estimation
estimate = c(para, em$mu0, em$Sigma0); SE = c(SE, NA, NA)
u = cbind(estimate, SE)
rownames(u) = c("phi","sigw","sigv","mu0","Sigma0"); u
Example 2.13
library(astsa)
set.seed(123)
num = 50
w = rnorm(num,0,.1)
x = cumsum(cumsum(w)) # states
y = x + rnorm(num,0,1) # observations
plot.ts(x, ylab="", lwd=2, ylim=c(-1,8))
lines(y, type='o', col='#808080')
# the model (specified by the parameters)
Phi = matrix(c(2,1,-1,0),2)
A = matrix(c(1,0),1)
mu0 = matrix(c(0,0),2)
Sigma0 = diag(c(1,1))
# innovations likelihood
Linn = function(para){
sigw = para[1]; sigv = para[2]; cQ = diag(c(sigw,0))
kf = Kfilter0(num,y,A,mu0,Sigma0,Phi,cQ,sigv)
return(kf$like)
}
# estimation
init.par = c(.1,1)
(est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE))
(SE = sqrt(diag(solve(est$hessian))))
# smooth and add to plot
sigw = est$par[1] # = 0.08 (se: 0.04)
cQ = diag(c(sigw,0))
sigv = est$par[2] # = 0.94 (se: 0.11)
ks = Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,sigv)
xsmoo = ts(ks$xs[1,1,]); psmoo = ts(ks$Ps[1,1,])
upp = xsmoo+2*sqrt(psmoo)
low = xsmoo-2*sqrt(psmoo)
lines(xsmoo, col=4, lty=2, lwd=3)
lines(upp, col=4, lty=2)
lines(low, col=4, lty=2)
# fit spline and add to plot
lines(smooth.spline(y), lty=1, col=2)
# add legends
legend("topleft", c("Observations","State"), pch=c(1,-1), lty=1, lwd=1:2, col=c('#808080',1))
legend("bottomright", c("Smoother", "GCV Spline"), lty=c(2,1), lwd=c(3,1), col=c(4,2) )
Example 2.15
library(astsa) data(WBC) num = length(WBC) A = array(ifelse(WBC>0, 1,0), dim=c(1,1,num)) # WBC = 0 when missing, otherwise, WBC > 0 # Run EM mu0 = 0; Sigma0 = .1; Phi = 1; cQ = .1; cR = .1 # initial parameter values (em = EM1(num,WBC,A,mu0,Sigma0,Phi,0,0,cQ,cR,0,100,.001)) # Run and graph smoother ks = Ksmooth1(num, WBC, A, em$mu0, em$Sigma0, em$Phi, 0, 0, sqrt(em$Q), sqrt(em$R), 0) Xs = ks$xs Ps = 3*sqrt(ks$Ps) plot(WBC, type="p", pch=20, ylim=c(1,5), xlab="day") lines(Xs) mss = 1 - as.vector(A) # mark missing data lines(mss, type="h", lwd=2) xx=c(time(WBC), rev(time(WBC))) # here and below puts conf intervals in gray yy=c(Xs-Ps, rev(Xs+Ps)) polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
Example 2.16
library(astsa)
data(jj)
num = length(jj)
A = cbind(1,1,0,0)
# Function to Calculate Likelihood
Linn=function(para){
Phi = diag(0,4); Phi[1,1] = para[1]
Phi[2,]=c(0,-1,-1,-1); Phi[3,]=c(0,1,0,0); Phi[4,]=c(0,0,1,0)
cQ1 = para[2]; cQ2 = para[3] # sqrt q11 and sqrt q22
cQ=diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
cR = para[4] # sqrt r11
kf = Kfilter0(num,jj,A,mu0,Sigma0,Phi,cQ,cR)
return(kf$like)
}
# Initial Parameters
mu0 = c(.7,0,0,0); Sigma0 = diag(.04,4)
init.par = c(1.03,.1,.1,.5) # Phi[1,1], the 2 Qs and R
# Estimation
est = optim(init.par,Linn,NULL,method="BFGS",hessian=TRUE,control=list(trace=1,REPORT=1))
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimate=est$par,SE)
rownames(u)=c("Phi11","sigw1","sigw2","sigv"); u
# Smooth
Phi = diag(0,4); Phi[1,1] = est$par[1]
Phi[2,]=c(0,-1,-1,-1); Phi[3,]=c(0,1,0,0); Phi[4,]=c(0,0,1,0)
cQ1 = est$par[2]; cQ2 = est$par[3]
cQ = diag(1,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
cR = est$par[4]
ks = Ksmooth0(num,jj,A,mu0,Sigma0,Phi,cQ,cR)
# Plot
Tsm = ts(ks$xs[1,,], start=1960, freq=4)
Ssm = ts(ks$xs[2,,], start=1960, freq=4)
p1 = 2*sqrt(ks$Ps[1,1,]); p2 = 2*sqrt(ks$Ps[2,2,])
dev.new(height=4)
par(mfrow=c(2,1), mar=c(3,3,1.5,1), mgp=c(1.6,.6,0),cex=.9)
plot(Tsm, main="Trend Component", ylab="Trend")
xx=c(time(jj), rev(time(jj)))
yy=c(Tsm-p1, rev(Tsm+p1))
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
plot(jj, main="Data and Trend+Season", ylab="J&J QE/Share")
xx = c(time(jj), rev(time(jj)))
yy=c( (Tsm+Ssm)-(p1+p2), rev((Tsm+Ssm)+(p1+p2)) )
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
# Forecast
n.ahead=12; y = ts(append(jj, rep(0,n.ahead)), start=1960, freq=4)
rmspe = rep(0,n.ahead); x00 = ks$xf[,,num]; P00 = ks$Pf[,,num]
Q=t(cQ)%*%cQ; R=t(cR)%*%(cR)
for (m in 1:n.ahead){
xp = Phi%*%x00; Pp = Phi%*%P00%*%t(Phi)+Q
sig = A%*%Pp%*%t(A)+R; K = Pp%*%t(A)%*%(1/sig)
x00 = xp; P00 = Pp-K%*%A%*%Pp
y[num+m] = A%*%xp; rmspe[m] = sqrt(sig)
}
dev.new(height=2.25)
par(mar=c(3,3,1.5,1), mgp=c(1.6,.6,0),cex=.9)
plot(y, type="o", main="", ylab="J&J QE/Share", ylim=c(5,30), xlim=c(1975,1984))
upp = ts(y[(num+1):(num+n.ahead)]+2*rmspe, start=1981, freq=4)
low = ts(y[(num+1):(num+n.ahead)]-2*rmspe, start=1981, freq=4)
lines(upp, lty=2); lines(low, lty=2); abline(v=1981, lty=3)
xx = c(time(low), rev(time(upp)))
yy = c(low, rev(upp))
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
[+] Chapter 3
Example 3.1
library(astsa) dev.new(height=4) par(mfrow=c(2,1), mar=c(2.5,2.8,.5,.5), mgp=c(1.6,.6,0)) plot(lynx/1000, ylab="lynx", type='o') plot(10*flu, col='#808080', ylab='flu mortality') Months=c( 'J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S', 'O', 'N', 'D') points(10*flu, pch=Months, cex=.7)
Example 3.2
library(astsa) par(mfrow=c(2,1)) plot(EEG) ar.fit = ar(EEG) # AIC best model p=35 innov = ar.fit$resid plot(innov, ylab = "Innovations") dev.new() par(mfrow=c(1,2)) acf(innov, na.action = na.omit) acf(innov^2, na.action = na.omit) dev.new() layout(matrix(c(1,2,4, 1,3,4), nc=2)) plot(sp500.gr, ylab="S&P500 Returns") acf(sp500.gr); acf(sp500.gr^2) plot(EXP6, ylab="Explosion")
Example 3.3
library(vcd) # load package - install it if you don't have it library(gamlss.data) # load package - ditto summary(goodfit(as.integer(discoveries), "nbinomial")) # GOF test summary(goodfit(as.integer(polio), "nbinomial")) # Graphic for the polio series - the discoveries graphic is similar layout(matrix(c(1,2, 1,3), nc=2)) plot(polio, type="h"); hist(polio); acf(polio)
Example 3.11
library(tsDyn) # load package - install it if you don't have it
library(astsa)
vignette("tsDyn") # for package details
lag1.plot(log10(lynx), 4, corr=FALSE)
(u = setar(log10(lynx), m=2, thDelay=1)) # fit model and view results
plot(u) # graphics - ?plot.setar for information
Example 3.12
library(fGarch) # load package - install it if you don't have it summary(fit <- garchFit(formula = ~arma(1,0)+aparch(1,1), data=sp500.gr)) v = volatility(fit, type = "sigma") par(mfrow=c(2,1)) plot(fit, which=c(1,2)) # use plot(fit) to see all options
Example 3.13
##-- Part 1 --## library(dyn) # load package - install it if you don't have it trend = time(polio) c1 = cos(2*pi*trend); s1 = sin(2*pi*trend) c2 = cos(2*pi*2*trend); s2 = sin(2*pi*2*trend) fit = dyn$glm(polio~trend+c1+s1+c2+s2+lag(polio,-1), family=quasipoisson, na.action=na.omit) summary(fit) # results plot(polio, type="h") lines(fitted(fit), col=4, lwd=2) ##-- Part 2 --## library(mgcv) # load package - install it if you don't have it fit2 = gam(polio ~ s(trend)+c1+s1+c2+s2, family="quasipoisson") summary(fit2) # results plot(polio, type="h") ffit2 = ts(fitted(fit2), start=1970, freq=12) lines(ffit2, col=4, lwd=2)
Example 3.14
library(nltsa) data(sp500.gr) bi.coh(sp500.gr) # Normalized BiSpectrum of the S&P500 returns dev.new() plot(sp500.gr) set.seed(666) # Normalized BiSpectrum of an AR(2) process (not in text) x.sim = arima.sim(list(order = c(2,0,0), ar = c(1,-.9)), n = 2048) dev.new() bi.coh(x.sim) dev.new() plot(x.sim)
Example 3.15
library(TSA) # load package - install it if you don't have it Keenan.test(log10(lynx)) Tsay.test(log10(lynx))
[+] Chapter 4
Figure 4.1
library(fGarch) # load package - install it if you don't have it set.seed(111) spec1 = garchSpec(model = list(omega=.3, alpha =.1, beta = .85)) spec2 = garchSpec(model = list(omega=.3, alpha =.1, beta = .899)) x1 = garchSim(spec1, n = 2000) x2 = garchSim(spec2, n = 2000) par(mfrow=c(2,1), mar=c(2,2,0,0)+.7, mgp=c(1.6,.6,0), oma=c(0,0,.75,0), cex=.9) plot.ts(x1) plot.ts(x2)
[+] Chapter 5
Figure 5.1
x0 = 0.3; nsmp = 100
x = matrix(0, nrow=nsmp, ncol=2)
u = runif(nsmp, min=0, max=1)
b = rbinom(nsmp, size=1, prob= 0.5)
x[1,1] = x0
# forward iteration
for (ismp in 2:nsmp){
x[ismp,1]= b[ismp]*x[ismp-1,1]*u[ismp]+(1-b[ismp])*((1-u[ismp])*x[ismp-1,1]+u[ismp])
}
# backward iteration
x[1,2]= x0
for (ismp in 2:nsmp){
x[ismp,2]= x0
for (j in 1:ismp-1){
x[ismp,2]= b[ismp-j]*x[ismp,2]*u[ismp-j]+(1-b[ismp-j])*((1-u[ismp-j])*x[ismp,2]+u[ismp-j])
}
}
par(mfrow=c(1,2), mar=c(3,2,2,1), mgp=c(1.6,.6,0), cex=.9)
plot(x[,1],type="l", ylim=c(0,1), ylab= "", main="Forward") # plot forward
plot(x[,2],type="l", ylim=c(0,1), ylab="", main="Backward") # plot backward
Example 5.45
library(logspline) # load package - install it if you don't have it
set.seed(90210); nsmp= 50; nbtraj= 10^4;
x = matrix(0, nrow=nsmp, ncol=nbtraj)
for (itraj in 1:nbtraj){
u = runif(nsmp,min=0,max=1); b = rbinom(nsmp, size=1, prob= 0.5)
x[1,itraj]= 0.01
for (ismp in 2:nsmp){ x[ismp,itraj]= b[ismp]*x[ismp-1,itraj]*u[ismp] +
(1-b[ismp])*((1-u[ismp])*x[ismp-1,itraj]+u[ismp])
}
}
z = seq(10^(-4),0.9999,length.out=100)
fz = 1/sqrt(z*(1-z))/pi
par(mfrow=c(2,2), mar=c(3,3,2,1))
plot(logspline(x[3,],lbound=0.0,ubound=1.0),lty=2,lwd=2,xlim=c(0,1),ylim=c(.5,3),main="x[3,]")
lines(z, fz, col="red")
plot(logspline(x[5,],lbound=0.0,ubound=1.0),lty=2,lwd=2,xlim=c(0,1),ylim=c(.5,3),main="x[5,]")
lines(z, fz, col="red")
plot(logspline(x[10,],lbound=0.0,ubound=1.0),lty=2,lwd=2,xlim=c(0,1),ylim=c(.5,3),main="x[10,]")
lines(z, fz, col="red")
plot(logspline(x[20,],lbound=0.0,ubound=1.0),lty=2,lwd=2,xlim=c(0,1),ylim=c(.5,3),main="x[20,]")
lines(z, fz, col="red")
Example 5.49
library(nltsa)
G = matrix(c(10,-5,-5,5), 2) # G is the target covariance
b = c(.1,3,100) # bG is the initial covariance
x0 = matrix(0,2)
niter = 10000
par(mfrow=c(3, 2))
for (i in 1:3){
out = metronorm(b[i]*G, G, niter, x0)
u = out$X[(niter-999):niter, 1]
plot.ts(u, ylab=expression(X[1]))
acf(u, 100, main="")
}
Example 5.52
set.seed(666)
y = rnorm(1000)
a0 = 3; b0 = 7; tau0 = 1; mu0 = 0; Nsim = 10000
ybar = mean(y); n = length(y)
an = a0 + n/2
tau = mu = rep(0, Nsim) # initialization arrays
tau[1] = rgamma(1, shape=a0, rate=b0) # sample an initial value
B = tau0 / (tau0 + n*tau[1])
mu[1] = rnorm(1, m=B*mu0+(1-B)*ybar, sd=1/sqrt(tau0+n*tau[1]))
for (t in 2:Nsim){
B = tau[t-1]/(tau0+n*tau[t-1])
mu[t] = rnorm(1, m=B*mu0+(1-B)*ybar, sd=1/sqrt(tau0+n*tau[t-1]))
ra1 = (1/2)*sum((y-mu[t])^2) + b0
tau[t] = rgamma(1,shape=an,rate=ra1)
}
par(mfrow=c(1,2)) # plot results
hist(mu[-1])
hist(tau[-1])
[+] Chapter 6
Example 6.16
set.seed(90210)
beta= 1.5
gamma= 0.25
nsmp = 20
nbtraj = 10^4
x = matrix(0,nrow=nsmp,ncol=nbtraj)
for (itraj in 1:nbtraj) { for (ismp in 2:nsmp)
{ x[ismp,itraj]= beta + gamma*( rpois(1,exp(x[ismp-1,itraj]))
- exp(x[ismp-1,itraj]))*exp(-x[ismp-1,itraj])
}
}
par(mfrow=c(2,2), mar=c(3,3,2,1))
xlow = beta-gamma; xup = beta+4*gamma
hist(x[3,], xlim=c(xlow,xup), freq=FALSE)
hist(x[5,], xlim=c(xlow,xup), freq=FALSE)
hist(x[10,], xlim=c(xlow,xup), freq=FALSE)
hist(x[20,], xlim=c(xlow,xup), freq=FALSE)
[+] Chapter 7
Figure 7.1
dev.new(height=5) par(mfcol=c(3,2), mar=c(3,3,1,1), mgp=c(1.6,.6,0), oma=c(0,0,1,0), cex=.8) set.seed(90210) phi = 0.95 x = arima.sim(list(ar=phi), n=100) y = cumsum(x)/(1:length(x)) cup = y + 1/sqrt(1:length(y))/(1-phi) clo = y - 1/sqrt(1:length(y))/(1-phi) plot(x, ylab=expression(X[~t])) mtext(expression(phi==+.95), side=3, line=.25) lmax = 20 ACF = as.vector(acf(x, lag.max = lmax, plot=FALSE)$acf) plot(0:lmax, ACF, type='h', xlab="Lag", ylim = c(-.1,1)) abline(h=0) trueacf = phi^(0:lmax) lines(0:lmax, trueacf, type="o", lty=2, col=2) plot(y, type='l', ylim=c(-5.0,10.0), xlab="n", ylab=expression(bar(X)[~n])) xx = c(1:length(x), length(x):1) yy = c(clo, rev(cup)) polygon(xx, yy, border=NA, col=gray(.6, alpha=.4)) abline(h=0, lty=2) phi = -.95 x = arima.sim(list(ar=phi), n=100) y = cumsum(x)/(1:length(x)) cup = y + 1/sqrt(1:length(y))/(1-phi) clo = y - 1/sqrt(1:length(y))/(1-phi) plot(x, ylab=expression(X[~t])) mtext(expression(phi==-.95), side=3, line=.25) lmax= 20 trueacf = phi^(0:lmax) ACF = as.vector(acf(x, lag.max = lmax, plot=FALSE)$acf) plot(0:lmax, ACF, type='h', xlab="Lag", ylim = c(-1,1)) abline(h=0) lines(0:lmax, trueacf, type="o", lty=2, col=2) plot(y, type="l", ylim=c(-.5,.5), xlab="n", ylab=expression(bar(X)[~n])); xx=c(1:length(x), length(x):1) yy=c(clo, rev(cup)) polygon(xx, yy, border=NA, col=gray(.6, alpha=.4)) abline(h=0, lty=2)
[+] Chapter 8
Example 8.30
library(nltsa) x = arima.sim(list(order=c(2,0,0), ar=c(1,-.9)), n=1000) # generate an ar2 u = arp.mcmc(x, porder=2, n.iter=1000, n.warmup=100) # sample from posterior # fun with the results library(MASS) # install it if it's not there already z = kde2d(u$phi[,1], u$phi[,2]) dev.new() plot(u$phi, xlab="phi1", ylab="phi2") contour(z, add=TRUE)
Example 8.31
library(BAYSTAR) # load package - install it if you don't have it x = log10(lynx) lagp1 = 1:2; lagp2 = 1:2 Iteration = 10000; Burnin = 2000 lynx.out = BAYSTAR(x,lagp1,lagp2,Iteration,Burnin,d0=2,step.thv=.5)
Example 8.32
library(astsa) library(bayesGARCH) # load package - install it if you don't have it library(IDPmisc) # load package - install it if you don't have it x = 100*diff(log(EuStockMarkets[,"CAC"])) y = window(x, start=1995, end=1998) plot(y); dev.new(); acf2(y); dev.new(); acf2(y^2) # not shown in text MCMC = bayesGARCH(as.vector(y)) plot(MCMC) # posterior statistics smpl = as.matrix(formSmpl(MCMC, l.bi = 2000)) ipairs(smpl)
[+] Chapter 9
Example 9.1
y = sqrt(EQcount) U = 2/sqrt(length(y)) lmax = 20 LAG = 1:lmax ACF = acf(y, plot=FALSE)$acf[-1] PACF = pacf(y, plot=FALSE)$acf par(mar=c(3,3,1,1), mgp=c(1.6,.6,0)) layout(matrix(c(1,2, 1,3), nc=2)) plot(EQcount, type='h') plot(LAG, ACF, type='h', ylim=c(-.2,.6)) abline(h=c(0,-U,U), lty=c(1,2,2), col=c(1,4,4)) plot(LAG, PACF, type='h', ylim=c(-.2,.6)) abline(h=c(0,-U,U), lty=c(1,2,2), col=c(1,4,4))
Example 9.2
library(mixtools) # load package - install it if you don't have it
# 3 normal mixture model
mixsp500 <- normalmixEM(sp500w.gr, lambda = c(0.4,0.4,0.2),
mu = c(0.00, 0.01,-0.01), sigma = c(0.01,0.01,0.01))
summary(mixsp500) # view results
# graphics
layout(matrix(1:4,2,byrow=TRUE), heights=rep(2:1,2))
par(mar=c(3,3,.5,.5), mgp=c(1.6,.6,0))
plot(ts(sp500w.gr,start=2003,freq=52), main="", ylab='S&P500 Weekly Returns')
hist(mixsp500$x, 25, prob=TRUE, main="", xlab="Data", xlim=c(-.25,.2), ylim=c(0,32))
for (i in 1:3) {mu = mixsp500$mu[i]; sig = mixsp500$sigma[i]
u = seq(-.2*i/3,.2*i/3, by=.001)
lines(u,dnorm(u, mean=mu, sd=sig), col=2*i+2)
}
# (P)ACF
U = 2/sqrt(length(sp500w.gr))
lmax = 30
LAG = 1:lmax
ACF = acf(sp500w.gr^2, lag.max= lmax, plot=FALSE)$acf[-1]
PACF = pacf(sp500w.gr^2, lag.max= lmax, plot=FALSE)$acf
plot(LAG, ACF, type='h', ylim=c(-.1,.3), main="")
abline(h=c(0,-U,U), lty=c(1,2,2), col=c(1,4,4))
plot(LAG, PACF, type='h', ylim=c(-.1,.3), main="")
abline(h=c(0,-U,U), lty=c(1,2,2), col=c(1,4,4))
# ### How the data were obtined ... uncomment and change dates if you want to update
# library(TTR) # install it if you don't have it
# Sys.setenv(tz='UTC')
# sp500 <- getYahooData("^GSPC",start=20021221,end=20120930,freq="daily")
# ep <- endpoints(sp500, on="weeks", k=1)
# sp500 = sp500[ep[-1]]
# sp500$sp500wr = diff(log(sp500$Close))
# sp500 = na.exclude(sp500)
# # for an xts plot
# plot(sp500$sp500wr, minor.ticks=FALSE, auto.grid=FALSE, main='')
# ### NOTE: sp500w.gr in 'nltsa.rda' is sp500$sp500wr
Example 9.6
library(nltsa) n = 100 u = NGMdata(n.obs=n, alpha=.5, beta=23, gamma=8, sig2v=1, sig2w=1, seed=90210) # sig2w will be 10 in Ch 12 y = u$y x = u$x # graphics dev.new(height=5) par(mfrow=c(2,2), mar=c(2.5,2.5,1,1), mgp=c(1.6,.6,0)) ts.plot(y, xlab="time", ylab=expression(y[~t])) ts.plot(x, xlab="time", ylab=expression(x[~t])) plot(x, y, xlab=expression(x[~t]), ylab=expression(y[~t])) plot.ts(x[-n], x[-1], ylab=expression(x[~t]), xlab=expression(x[~t-1]), cex=.8)
Example 9.11
library(depmixS4)
library(nltsa)
## estimation ##
model <- depmix(EQcount ~1, nstates=2, data=data.frame(EQcount), family=poisson())
set.seed(90210)
(fm <- fit(model)) # fm contains results from estimation
summary(fm)
#-- get parameters --#
# note- can't control the state names in depmix, so here we
# maintain which is state 1 [min lam] and which is state 2 [max lam]
u = as.vector(getpars(fm))
if (u[7] <= u[8]) { para.mle = c(u[3:6], exp(u[7]), exp(u[8]))
} else { para.mle = c(u[6:3], exp(u[8]), exp(u[7]))
}
mtrans = matrix(para.mle[1:4], byrow=TRUE, nrow=2)
lams = para.mle[5:6]
pi1 = mtrans[2,1]/(2 - mtrans[1,1] - mtrans[2,2])
pi2 = 1 - pi1
## graphics ##
dev.new(height=4)
layout(matrix(c(1,2,1,3), 2))
par(mar=c(3,3,1,1), mgp=c(1.6,.6,0))
# data and states
plot(EQcount, main="", ylab='EQcount', type='h', col='#808080')
text(EQcount, col=posterior(fm)[,1], labels=posterior(fm)[,1], cex=.9)
# prob of state 2
plot(ts(posterior(fm)[,3], start=1900), ylab=expression(hat(pi)[~t]*' (2 | data)'))
abline(h=.5, lty=2)
# histogram
hist(EQcount, breaks=30, prob=TRUE, main="")
xvals = seq(1,45)
u1 = pi1*dpois(xvals, lams[1])
u2 = pi2*dpois(xvals, lams[2])
lines(xvals, u1, col=3)
lines(xvals, u2, col=4)
### Bootstap ####
# function to generate data
pois.HMM.generate_sample = function(n,m,lambda ,Mtrans,StatDist=NULL){
# n = data length, m = # of states, Mtrans = transition matrix, StatDist = stationary distn
if(is.null(StatDist))StatDist =solve(t(diag(m)-Mtrans +1),rep(1,m))
mvect = 1:m
state = numeric(n)
state [1] = sample(mvect ,1,prob=StatDist)
for (i in 2:n)
state[i]=sample(mvect ,1,prob=Mtrans[state[i-1] ,])
y = rpois(n,lambda=lambda[state ])
list(y= y, state= state)
}
# start it up
set.seed(10101101)
nboot = 100
nobs = length(EQcount)
para.star = matrix(NA, nrow=nboot, ncol = 6)
for (j in 1:nboot){
x.star = pois.HMM.generate_sample(n=nobs, m=2, lambda=lams, Mtrans=mtrans)$y
model <- depmix(x.star ~1, nstates=2, data=data.frame(x.star), family=poisson())
u = as.vector(getpars(fit(model, verbose=0)))
# make sure state 1 is the one with the smaller intensity parameter
if (u[7] <= u[8]) { para.star[j,] = c(u[3:6], exp(u[7]), exp(u[8])) }
else { para.star[j,] = c(u[6:3], exp(u[8]), exp(u[7])) }
}
# bootstrapped stnd errors
SE = sqrt(apply(para.star,2,var)+ (apply(para.star,2,mean)-para.mle)^2)[c(1,4:6)]
names(SE)=c('seM11/M12', 'seM21/M22', 'seLam1', 'seLam2')
SE
Example 9.12
library(depmixS4)
library(nltsa)
x = ts(sp500w.gr, start=2003, freq=52) # cleans up the data a bit
# fit 3-state model
mod3 <- depmix(x~1, nstates=3, data=data.frame(x))
set.seed(2)
(fm3 <- fit(mod3))
# graphics
dev.new(height=3)
par(mar=c(3,3,1,1), mgp=c(1.6,.6,0))
plot(x, main="", ylab='S&P500 Weekly Returns', col='#9f9f9f', ylim=c(-.11,.11))
text(x, col=4-posterior(fm3)[,1], labels=4-posterior(fm3)[,1], cex=.9)
# NOTE: For display purposes, we reversed the names of states 1 and 3,
# the graphic looks better this way.
# the MLEs:
para.mle = as.vector(getpars(fm3)[-(1:3)])
# for display (states 1 and 3 names switched)
permu = matrix(c(0,0,1,0,1,0,1,0,0), 3,3)
(mtrans.mle = permu%*%round(t(matrix(para.mle[1:9],3,3)),3)%*%permu)
(norms.mle = permu%*%round(matrix(para.mle[10:15], 3,2),3))
# bootstrap
set.seed(90210)
n.obs = length(x)
n.boot = 200 # lower this if it takes too long
para.star = matrix(NA, nrow=n.boot, ncol = 15)
respst <- para.mle[10:15]
trst <- para.mle[1:9]
for ( nb in 1:n.boot){
mod <- simulate(mod3)
x.star = as.vector(mod@response[[1]][[1]]@y)
dfx = data.frame(x.star)
mod.star <- depmix(x.star~1, data=dfx, respst=respst, trst=trst, nst=3)
fm.star = fit(mod.star, verbose=0)
para.star[nb,] = as.vector(getpars(fm.star)[-(1:3)])
}
# bootstrap stnd errors
SE = sqrt(apply(para.star,2,var)+ (apply(para.star,2,mean)-para.mle)^2)
# for display
(SE.mtrans.mle = permu%*%round(t(matrix(SE[1:9],3,3)),3)%*%permu)
(SE.norms.mle = permu%*%round(matrix(SE[10:15], 3,2),3))
Example 9.13
library(astsa) data(flu) library(MSwM) # fit model on diffed flu dflu = as.vector(diff(flu)) model = lm(dflu~1) mod = msmFit(model, k=2, p=2, sw=rep(TRUE,4)) # 2 regimes, AR(2)s summary(mod) plotProb(mod) # smoother/filter probs for state 2 sp = mod@Fit@smoProb[,2] fp = mod@Fit@filtProb[,2] # graphic x = window(diff(flu), start= c(1968,3)) regime = rep(2, length(x)) regime[sp <.5] = 1 # for graphic, get regime numbers 1 or 2 dev.new(height=3) layout(matrix(1:2,2,2), heights=c(2.2,1)) par(mar=c(2, 3,.2,.5), mgp=c(1.6,.6,0), cex=.8) plot(x, col='#808080', ylab= expression(nabla~flu[~t]), xlab='') abline(v=1968:1978,lty=3, col=8); abline(h=seq(-.4,.4,.2), lty=3, col=8) text(x, col=regime, labels=regime, cex=.9) sp = ts(sp, start=tsp(x)[1], freq=12) fp = ts(fp, start=tsp(x)[1]+1/12, freq=12) plot(sp, col=1, ylab= expression(hat(pi)[~t]*' (2 | n)'), axes=FALSE, xlab='') axis(1); axis(2, c(0,.5,1)); box() abline(v=1968:1978,lty=3, col=8); abline(h=c(.25,.5,.75), lty=3, col=8) lines(fp, type='h', col=4)
[+] Chapter 10
Example 10.1
## Fig 10.1 ##
layout(matrix(c(1:12), 3,4, byrow = TRUE))
set.seed(3142)
# parameters of the Gaussian mixture
mu0= -1; sigma0= 2.5; mu1= 1.5; sigma1= 0.3; alpha= 0.9
df= 4 # df of the t-distribution
nbsmp= 10^3 # number of samples
indx = sample(1:nbsmp, 150) # indeces for displaying weights
par(mar=c(1,1,.5,0)+1, mgp=c(1.6,.6,0), oma=c(0,1.5,.5,0))
cnt = 0
for (scalet in c(.5,5,15)) {
x = seq(-10,10,length.out=300)
lx = length(x);
y = alpha*dnorm(x,mu0*matrix(1,nrow=1,ncol=lx),
+ sigma0*matrix(1,nrow=1,ncol=lx))+ (1-alpha)*dnorm(x,mu1*matrix(1,nrow=1,ncol=lx),
+ sigma1*matrix(1,nrow=1,ncol=lx))
z = dt(x/scalet,df*matrix(1,nrow=1,ncol=lx))/scalet
plot(x,y,type='l',col=1, ylim=c(0, .2),lwd=2,ylab="",xlab="")
lines(x,z,lty=2,col=4,lwd=2)
mtext(paste("scale =", scalet), side=2, line=2)
if (cnt < 1) mtext("Target/Proposal", line=.5)
# Computation of the importance function
w = y/z
plot(x,w,type='l',lwd=2, ylab="", xlab="")
if (cnt < 1) mtext("Importance Function", line=.5)
# Sample the proposal distribution
randsmp = scalet*rt(nbsmp,df);
wrnum = alpha*dnorm(randsmp,mu0*matrix(1,nrow=1,ncol=nbsmp),
+ sigma0*matrix(1,nrow=1,ncol=nbsmp))+ (1-alpha)*dnorm(randsmp,mu1*matrix(1,nrow=1,ncol=nbsmp),
+ sigma1*matrix(1,nrow=1,ncol=nbsmp))
wrden = dt(randsmp/scalet,df*matrix(1,nrow=1,ncol=nbsmp))/scalet
# Computation of the importance weights
wr = wrnum / wrden
plot(randsmp[indx], wr[indx], type='h', xlim=c(-10.0,10.0), ylab="", xlab="")
if (cnt < 1) mtext("Weights", line=.5)
# computation of the estimator
plot(cumsum(randsmp * wr)/cumsum(wr), type= 'l', lwd=2, ylab= "", ylim=c(-1.5,0), xlab= "")
if (cnt < 1) mtext("Estimator", line=.5)
cnt=cnt+1
}
Example 10.8
##################
##-- Fig 10.2 --##
##################
library(ggplot2)
library(nltsa)
library(reshape2)
set.seed(19091983)
npart <- 5000
arch <- ARCH(alpha0=1, alpha1=0.99,sigmaV=1)
simulated.data <- simulate.data(arch, x1=0, 20)
true_state <- simulated.data$x
T <- length(true_state)
all.particles <- data.frame()
for (sigmaV in c(1,10)) {
arch$sigmaV <- sigmaV
# Prior kernel
b.prior <- siskernel(arch, true_state, npart,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Optimal kernel
b.optimal <- siskernel(arch, true_state, npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
p.prior <- melt(t(b.prior$particles), varnames=c("Index","Iteration"), value.name="Particles")
p.optimal <- melt(t(b.optimal$particles), varnames=c("Index","Iteration"), value.name="Particles")
p.prior$Kernel <- as.factor('Prior')
p.prior$sigmaV <- sigmaV
p.optimal$Kernel <- as.factor('Optimal')
p.optimal$sigmaV <- sigmaV
all.particles <- rbind(all.particles, p.prior, p.optimal)
}
dev.new(height=4)
all.particles$Iteration <- as.factor(all.particles$Iteration)
observations <- data.frame(Iteration=1:length(true_state),State=true_state)
qplot(Iteration,Particles,geom="violin",facets=sigmaV~Kernel,data=all.particles) +
theme_bw() +
facet_grid(sigmaV~Kernel, labeller="label_both") +
geom_line(data=observations,aes(Iteration,State),color="red",alpha=.5) +
geom_point(data=observations,aes(Iteration,State),color="red",shape=4) +
coord_cartesian(ylim = c(-21.2, 21.2)) + scale_x_discrete(breaks=seq(0,20,by=2))
##################
##-- Fig 10.3 --##
##################
library(ggplot2)
library(nltsa)
library(reshape2)
set.seed(19091983)
ARCH.compare.mean <- function(arch, y, nrep=100,npart=5000) {
T <- length(y)
b.prior <- sis(arch, y, N=npart)
b.optimal <- siskernel(arch, y, N=npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf)
estim.prior <- matrix(nrow=nrep,ncol=T)
estim.optimal <- matrix(nrow=nrep,ncol=T)
results <- data.frame()
begin <- Sys.time()
for (k in seq(1,nrep)) {
# Prior kernel
b.prior <- siskernel(arch, y, npart,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
estim.prior[k,] <- rowSums(b.prior$particles * b.prior$weights)
# Optimal kernel
b.optimal <- siskernel(arch, y, npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
estim.optimal[k,] <- rowSums(b.optimal$particles * b.optimal$weights)
# Estimate remaining runtime
now <- Sys.time()
will.end <- now+(nrep-k)*(now - begin)/k
print(sprintf("Iteration %d/%d, expected end %s",k,nrep,format(will.end)))
}
center=FALSE
Estimates <- c(as.vector(scale(estim.prior,center=center,scale=FALSE)),
as.vector(scale(estim.optimal,center=center,scale=FALSE)))
Time <- as.factor(rep(b.prior$t, each=nrep, times=2))
Kernel <- factor(rep(c("Prior","Optimal"),each=nrep*T),ordered=TRUE)
data.frame(Time,Estimates,Kernel)
}
arch <- ARCH(alpha0=1, alpha1=0.99, sigmaV=1)
simulated.data <- simulate.data(arch, x1=0, 20)
true_state <- simulated.data$x
T <- length(true_state)
arch$sigmaV <- 10
large.sigma <- cbind(ARCH.compare.mean(y=true_state, arch),sigmaV=arch$sigmaV)
arch$sigmaV <- 1
small.sigma <- cbind(ARCH.compare.mean(y=true_state, arch),sigmaV=arch$sigmaV)
all.sigma <- rbind(large.sigma,small.sigma)
# graphic
dev.new(height=4)
observations <- data.frame(Time=1:length(true_state),State=true_state)
qplot(Time,Estimates,geom="violin",facets=sigmaV~Kernel,data=all.sigma) +
theme_bw() +
facet_grid(sigmaV~Kernel, labeller="label_both") +
geom_line(data=observations,aes(Time,State),color="red",alpha=.5) +
geom_point(data=observations,aes(Time,State),color="red",shape=4) +
scale_x_discrete(breaks=seq(0,20,by=2))
Example 10.9
##################
##-- Fig 10.4 --##
##################
library(ggplot2)
library(nltsa)
library(ineq)
library(reshape2)
set.seed(19091983)
arch <- ARCH(alpha0=1,alpha1=0.99,sigmaV=1)
npart <- 5000
T <- 100
simulated.data <- simulate.data(arch, x1=0, T)
y <- simulated.data$y
# Prior kernel
b.prior <- siskernel(arch, y, npart,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Optimal kernel
b.optimal <- siskernel(arch, y, npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Plot Lorenz curve of the importance weights
IterationOfInterest <- c(5,10,25,50)
results <- data.frame()
for (i in IterationOfInterest) {
x <- Lc(b.prior$weights[i,])
y <- data.frame(p=x$p, L=x$L, Iteration=i, Kernel=as.factor('Prior'))
results <- rbind(results, y)
}
for (i in IterationOfInterest) {
x <- Lc(b.optimal$weights[i,])
y = data.frame(p=x$p, L=x$L, Iteration=i, Kernel=as.factor('Optimal'))
results <- rbind(results, y)
}
# graphic
dev.new(height=3)
ggplot(aes(100*p, L), data=results) +
facet_grid(Kernel~Iteration, labeller="label_both") +
ylab("Total weight") +
xlab("Percentage of particles") +
coord_cartesian(xlim = c(-10, 110)) +
theme_bw() +
geom_path()
##################
##-- Fig 10.5 --##
##################
set.seed(19091983)
arch <- ARCH(alpha0=1, alpha1=0.99, sigmaV=1)
npart <- 5000
T <- 100
simulated.data <- simulate.data(arch, x1=0, T)
y <- simulated.data$y
# Prior kernel
b.prior <- siskernel(arch, y, npart,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Optimal kernel
b.optimal <- siskernel(arch, y, npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
## Display histograms of weights on log scale for optimal kernel ##
# extract data
w.prior <- melt(t(b.prior$weights), varnames=c("Particle","Time"), value.name="Weights")
w.prior$Kernel <- factor('Prior')
w.optimal <- melt(t(b.optimal$weights), varnames=c("Particle","Time"), value.name="Weights")
w.optimal$Kernel <- factor('Optimal')
w.both <- rbind(w.prior, w.optimal)
# graphic
clr = 'lightblue'
brks = seq(-10,-1,by=.05)
dev.new(height=3.5)
par(mfrow=c(4,1), mar=c(1.5,3,1,.5), mgp=c(1.6,.6,0))
for ( i in c(5,10,25,50) ){
datax = subset(w.both, Time %in% i & Kernel %in% c('Optimal'))
u = hist(log10(datax[,'Weights']), breaks=brks, col=clr, xlim=c(-10,-1.5), xlab="", main="",
panel.first=grid(ny=6, col=8))
legend('topleft', paste(' Iteration:', i), bty='n', cex=1.25 )
}
Figure10.6
library(ggplot2)
library(nltsa)
library(reshape2)
theme_set(theme_bw())
set.seed(19091983)
arch <- ARCH(alpha0=1, alpha1=0.99,sigmaV=1)
npart <- 5000
T <- 100
simulated.data <- simulate.data(arch, x1=0, T)
y <- simulated.data$y
# Prior kernel
b.prior <- siskernel(arch, y, npart,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Optimal kernel
b.optimal <- siskernel(arch, y, npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
timeOfInterest <- 1:T
results <- data.frame()
for (t in timeOfInterest) {
x <- ESS(b.prior$weights[t,])
y <- data.frame(ESS=x, Timestep=t, Kernel=as.factor('Prior'))
results <- rbind(results, y)
}
for (t in timeOfInterest) {
x <- ESS(b.optimal$weights[t,])
y = data.frame(ESS=x, Timestep=t, Kernel=as.factor('Optimal'))
results <- rbind(results, y)
}
# graphic
dev.new(height=2.5)
par(mar=c(3,3,1,1), mgp=c(1.6,.6,0), cex=.9)
yP = results[1:100,]
yO = results[101:200,]
plot.ts(yO$ESS, xlab='Time', ylab="Effective Sample Size", log='y', ylim=c(1,5000),
panel.first=grid(col=8))
lines(yP$ESS, col=4, lwd=2)
text(x=40,y=1200, "Optimal Kernel", cex=.9)
text(x=16,y=50, "Prior Kernel", col=4, cex=.9)
Example 10.13
## Figure 10.7 ##
library(ggplot2)
library(nltsa)
library(reshape2)
library(ineq)
set.seed(19091983)
arch <- ARCH(alpha0=1,alpha1=0.99,sigmaV=1)
npart <- 5000
T <- 100
simulated.data <- simulate.data(arch, x1=0, T)
y <- simulated.data$y
# Prior kernel
b.prior <- sisr(arch, y, npart,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Optimal kernel
b.optimal <- sisr(arch, y, npart,
proposal.rnd=ARCH.optimal.rnd,
proposal.logpdf=ARCH.optimal.logpdf,
initial.proposal.rnd = ARCH.optimal.initial.rnd,
initial.proposal.logpdf = ARCH.optimal.initial.logpdf)
# Display histograms of weifghts on log scale for optimal kernel
w.prior <- melt(t(b.prior$weights), varnames=c("Particle","Time"), value.name="Weights")
w.prior$Kernel <- factor('Prior')
w.optimal <- melt(t(b.optimal$weights), varnames=c("Particle","Time"), value.name="Weights")
w.optimal$Kernel <- factor('Optimal')
w.both <- rbind(w.prior, w.optimal)
# graphic
clr = 'lightblue'
brks = seq(-5,-3.5,by=.01)
dev.new(height=3.5)
par(mfrow=c(4,1), mar=c(1.5,3,1,3), mgp=c(1.6,.6,0))
for ( i in c(5,10,25,50) ){
datax = subset(w.both, Time %in% i & Kernel %in% c('Optimal'))
u = hist(log10(datax[,'Weights']), breaks=brks, col=clr, xlim=c(-5,-3.5), xlab="", main="",
panel.first = grid(ny=6, col=8))
legend('topleft', paste(' Iteration:', i), bty='n', cex=1.25 )
}
Figure 10.8
# parameters of the Gaussian mixture
mu0= -1; sigma0= 2.5; mu1= 1.5; sigma1= 0.3; alpha= 0.9;
xmin=-10; xmax=10;
# parameters of the t-distribution
df= 4;
x= seq(xmin,xmax,length.out=300)
lx= length(x); deltax= (xmax-xmin)/lx
y= alpha*dnorm(x,mu0*matrix(1,nrow=1,ncol=lx),
sigma0*matrix(1,nrow=1,ncol=lx))+ (1-alpha)*dnorm(x,mu1*matrix(1,nrow=1,ncol=lx),
sigma1*matrix(1,nrow=1,ncol=lx))
scalet= seq(0.1,20,length.out=100)
lscalet= length(scalet)
variance= rep(0,lscalet)
for (iscale in 1:lscalet){
scalecur= scalet[iscale]
# Computation of the importance function
z = dt(x/scalecur,df*matrix(1,nrow=1,ncol=lx))/scalecur
w = y/z
# Computation of the variance
variance[iscale]= deltax*mean(z*w^2*((x-alpha*mu0-(1-alpha)*mu1)^2))
}
# graphic
dev.new(height=2.25)
par(mar=c(3,3,1,1), mgp=c(1.6,.6,0), cex=.95)
plot(scalet,log(variance,base=10), type='l', xlab= "Scale", ylab= expression(log[10]~Variance),
panel.first= grid(col=8))
[+] Chapter 11
Example 11.2
library(ggplot2)
library(plyr)
library(reshape2)
library(testthat)
library(nltsa)
set.seed(90210)
T <- 40
n.particles <- 50
nlss <- NoisyAR(phi=0.9, sigmaW=.1, sigmaV=.2)
data <- simulate.data(nlss, x1=0, T=T)
truex.df <- melt(data$x, varnames=c("Time",""))
# Call the poor man's smoother
p <- sisr(nlss, y=data$y, N=n.particles)
smoothed.poor <- poorman.smoother(p)
# See that filtering and smoothing are equivalent for the last time step
expect_equal(smoothed.poor$weights,p$weights[T,])
expect_equal(smoothed.poor$particles[T,,],p$particles[T,,])
# Recover unweighted sample by resampling
indices.poor <- ResidualMultinomialR(smoothed.poor$weights)
resampled.poor <- smoothed.poor$particles[,indices.poor,,drop=FALSE]
# Plot trajectories
traj.poor.df <- data.frame(melt(resampled.poor, varnames=c("Time","Index","Dimension")),
Algorithm="Poorman")
dev.new(height=2.25)
ggplot(data=traj.poor.df, aes(x=factor(Time),y=value)) +
geom_line(aes(group=as.factor(Index), colour=as.factor(Index))) +
facet_grid(Algorithm~.) + theme_bw() +
scale_fill_grey() + scale_colour_grey() +
ylab("smoothed particles") + xlab("time") +
theme(legend.position='none') +
scale_x_discrete(breaks=seq(0,40,by=5))
Example 11.7
library(nltsa)
library(ggplot2)
library(plyr)
library(reshape2)
library(astsa)
#######################
##### Fig 11.2 ########
#######################
set.seed(19091983)
T <- 200
n.particles <- 100
nlss <- NoisyAR(phi=0.8, sigmaW=.1, sigmaV=.2)
data <- simulate.data(nlss,x1=0,T=T)
truex.df <- melt(data$x,varnames=c("Time",""))
#======= Study the distribution of the particles ===========
# Call the poor man's smoother
p <- sisr(nlss, y=data$y, N=n.particles)
smoothed.poor <- poorman.smoother(p)
smoothed.ffbsm <- FFBSm(nlss,p)
smoothed.ffbsi <- FFBSi(nlss,p)
smoothed.ffbsilin <- FFBSi.linear(nlss,p,sigma.plus=1/(nlss$sigmaW*sqrt(2*pi)))
# PGAS with n.particles MCMC steps and only 10 particles at each MCMC iteration
smoothed.pgas <- PGAS(nlss=nlss, y=data$y, n.particles=10, n.mcmc=n.particles)
# Recover unweighted sample by resampling
indices.poor <- ResidualMultinomialR(smoothed.poor$weights)
resampled.poor <- smoothed.poor$particles[,indices.poor,,drop=FALSE]
resampled.ffbsm <- array(dim=dim(smoothed.ffbsm$particles))
for (k in 1:T) {
indices.ffbsm <- ResidualMultinomialR(smoothed.ffbsm$weights[k,])
resampled.ffbsm[k,,] <- smoothed.ffbsm$particles[k,indices.ffbsm,,drop=FALSE]
}
resampled.ffbsi <- smoothed.ffbsi$particles
resampled.ffbsilin <- smoothed.ffbsilin$particles
resampled.pgas <- smoothed.pgas$particles
# Plot remaining trajectories (after a final resampling to ensure uniform weights)
traj.poor.df <- data.frame(melt(resampled.poor, varnames=c("Time","Index","Dimension")),
Algorithm="Poorman")
traj.ffbsm.df <- data.frame(melt(resampled.ffbsm, varnames=c("Time","Index","Dimension")),
Algorithm="FFBSm")
traj.ffbsi.df <- data.frame(melt(resampled.ffbsi, varnames=c("Time","Index","Dimension")),
Algorithm="FFBSi")
traj.ffbsilin.df <- data.frame(melt(resampled.ffbsilin, varnames=c("Time","Index","Dimension")),
Algorithm="FFBSi Linear")
traj.pgas.df <- data.frame(melt(resampled.pgas, varnames=c("Time","Index","Dimension")),
Algorithm="PGAS")
traj.df <- rbind(traj.poor.df,
traj.ffbsm.df,
traj.ffbsi.df,
traj.ffbsilin.df,
traj.pgas.df)
# Compute the quantiles
q <- ddply(traj.df, ~Algorithm+Time,
function (x) quantile(x$value, probs = c(.05, .5, .95), names=FALSE))
names(q)[3] <- "q05"
names(q)[4] <- "q50"
names(q)[5] <- "q95"
ksmooth <- Ksmooth0.NoisyAR(nlss,data$y)
xs <- ksmooth$xs
xf <- ksmooth$xf
xs.df <- melt(xs,varnames=c("","","Time"))
xf.df <- melt(xf,varnames=c("","","Time"))
ggplot(data=q,aes(x=Time,y=q50))+
ylab("5%-95% quantiles of smoothed particles") +
xlab("time") + geom_line(linetype='dashed', alpha=.5) +
geom_ribbon(aes(ymin=q05,ymax=q95), alpha=.5) +
facet_grid(Algorithm~.) + theme_bw() +
scale_color_grey() +
theme(legend.position="none") +
geom_line(data=xs.df,aes(group=1,y=value),linetype='solid')
#######################
##### Fig 11.3 ########
#######################
set.seed(19091983)
T <- 200
n.particles <- 100
nrep <- 20
nlss <- NoisyAR(phi=0.8,sigmaW=.1,sigmaV=.2)
data <- simulate.data(nlss,x1=0,T=T)
truex.df <- melt(data$x,varnames=c("Time",""))
ksmooth <- Ksmooth0.NoisyAR(nlss,data$y)
xs <- ksmooth$xs
xf <- ksmooth$xf
xs.df <- melt(xs,varnames=c("","","Time"))
xf.df <- melt(xf,varnames=c("","","Time"))
# Compute nrep posterior mean estimates: a bit long to run
estimates.poor <- array(data=NA_real_, dim=c(nrep,T) )
estimates.ffbsm <- array(data=NA_real_, dim=c(nrep,T))
estimates.ffbsi <- array(data=NA_real_, dim=c(nrep,T))
estimates.ffbsilin <- array(data=NA_real_, dim=c(nrep,T))
estimates.pgas <- array(data=NA_real_, dim=c(nrep,T))
pr <- progress_text()
pr$init(nrep*T)
for (r in 1:nrep) {
bstrp <- sisr(nlss,y=data$y,N=n.particles)
poor <- poorman.smoother(bstrp)
ffbsm <- FFBSm(nlss,bstrp)
ffbsi <- FFBSi(nlss,bstrp)
ffbsilin <- FFBSi.linear(nlss,bstrp,sigma.plus=1/(nlss$sigmaW*sqrt(2*pi)))
pgas <- PGAS(nlss=nlss, y=data$y, n.particles=15, n.mcmc=n.particles)
for (k in T:1) {
estimates.poor[r,k] <- sum(poor$particles[k,,]*poor$weights/sum(poor$weights))
estimates.ffbsm[r,k] <- sum(ffbsm$particles[k,,]*ffbsm$weights[k,]/sum(ffbsm$weights[k,]))
estimates.ffbsi[r,k] <- mean(ffbsi$particles[k,,])
estimates.ffbsilin[r,k] <- mean(ffbsilin$particles[k,,])
estimates.pgas[r,k] <- mean(pgas$particles[k,,])
pr$step()
}
}
# Plot the independent estimated smoothed trajectories
est.poor.df <- data.frame(melt(estimates.poor,varnames=c("Repetition","Time")),
Algorithm="Poorman")
est.ffbsm.df <- data.frame(melt(estimates.ffbsm,varnames=c("Repetition","Time")),
Algorithm="FFBSm")
est.ffbsi.df <- data.frame(melt(estimates.ffbsi,varnames=c("Repetition","Time")),
Algorithm="FFBSi")
est.ffbsilin.df <- data.frame(melt(estimates.ffbsilin,varnames=c("Repetition","Time")),
Algorithm="FFBSi Linear")
est.pgas.df <- data.frame(melt(estimates.ffbsilin,varnames=c("Repetition","Time")),
Algorithm="PGAS")
est.df <- rbind(est.poor.df,
est.ffbsm.df,
est.ffbsi.df,
est.ffbsilin.df,
est.pgas.df)
# Compute the quantiles
q <- ddply(est.df, ~Algorithm+Time, function (x) quantile(x$value,
probs = c(.05, .5, .95), names=FALSE))
names(q)[3] <- "q05"
names(q)[4] <- "q50"
names(q)[5] <- "q95"
ksmooth <- Ksmooth0.NoisyAR(nlss,data$y)
xs <- ksmooth$xs
xf <- ksmooth$xf
xs.df <- melt(xs,varnames=c("","","Time"))
xf.df <- melt(xf,varnames=c("","","Time"))
# save(list = ls(all=TRUE),file = "smooth.Rdata") # if you want to save the output
ggplot(data=q,aes(x=Time,y=q50))+
ylab("5%-95% quantiles of smoothed mean estimates") + xlab("time") +
geom_line(linetype='dashed', alpha=.5) +
geom_ribbon(aes(ymin=q05,ymax=q95), alpha=.5) +
facet_grid(Algorithm~.) + theme_bw() +
scale_color_grey() +
theme(legend.position="none") +
geom_line(data=xs.df,aes(group=1,y=value),linetype='solid')
[+] Chapter 12
Example 12.5
library(nltsa)
library(reshape2)
library(plyr)
# n.rep = 10 and n.em = 25 or 30 is ok and is much quicker
n.rep = 25
n.em = 50
n.particles = function(i) min(50+i*2, 500)
compute.time = TRUE # do you want to know how long this will take you (FALSE=no)
T = 100
alpha = .5
beta = 25
gamma = 8
sig2w = 10
sig2v = 1
mcmseed = 90210
nlss = NGMmodel(alpha, beta, gamma, sig2w, sig2v)
data = NGMdata(n.obs=T, alpha, beta, gamma, sig2w=10, sig2v=1, seed=mcmseed) # the default
obs = as.matrix(data$y) # ParticleEM needs a vector
######
if (compute.time)
{ singleTime <- system.time(ParticleEM(initial.nlss=nlss, y=obs,
n.particles=n.particles(n.em), n.em=1))[3]
cat('/n')
print(sprintf('One single EM iteration with %d particles takes %.2f seconds',
n.particles(n.em), singleTime))
hour <- floor(n.em*singleTime/3600)
min <- floor(n.em*singleTime/60 - hour*60)
sec <- n.em*singleTime - min*60 - hour*3600
print(sprintf('Estimated runtime for %d EM iterations: %dh %dm %.2fs',
n.em, hour, min, sec))
}
######
alltraces = array(NA, dim=c(n.em+1, 5, n.rep))
ptm <- proc.time()
set.seed(mcmseed)
for (irep in 1:n.rep){
alpha = runif(n=1, min=0.5/2, max=0.5*1.5)
beta = runif(n=1, min=25/2, max=25*1.5)
gamma = runif(n=1, min=8/2, max=8*1.5)
sigmaW = runif(n=1, min=sqrt(10)/2, max=sqrt(10)*1.5)
sigmaV = runif(n=1, min=1/2, max=1*1.5)
initial.nlss = NGMmodel(alpha, beta, gamma, sig2w=sigmaW^2, sig2v=sigmaV^2)
alltraces[,,irep]= ParticleEM(initial.nlss=initial.nlss, y = obs,
n.particles=n.particles, n.em=n.em)$trace
cat("\n")
}
(time2run = proc.time() - ptm)
# graphics
truep = c(.5,25,8,sqrt(10),1)
names = c(expression(alpha), expression(beta), expression(gamma),
expression(sigma[w]), expression(sigma[v])) # ParticleEM returns stnd. devs
dev.new(height=3)
par(mfcol=c(1,5), mar=c(3,1,1,1), mgp=c(1.6,.6,0), oma=c(0,1.5,1,0))
for (i in 1:5){
ts.plot(alltraces[,i,], ylab='', xlab='iteration')
mtext(names[i], line=.25, side=3, cex=1.1)
if (i == 1) mtext('trace', side=2, line=1.5, cex=.9)
abline(h=truep[i], col=2, lwd=1)
}
Example 12.6
library(nltsa)
library(reshape2)
library(plyr)
# n.rep = 10, n.em = 500, and n.particles = 5 is ok if this takes too long
n.rep = 30
n.em = 1000
n.particles = 15
T = 100
alpha = .5
beta = 25
gamma = 8
sig2w = 10
sig2v = 1
mcmseed = 90210
nlss = NGMmodel(alpha, beta, gamma, sig2w, sig2v)
data = NGMdata(n.obs=T, alpha, beta, gamma, sig2w=10, sig2v=1, seed=mcmseed)
obs = as.matrix(data$y)
alltracesAve = array(NA, dim=c(n.em+1, 5, n.rep))
ptm <- proc.time()
for (irep in 1:n.rep){
# Draw starting values of the algorithm
alpha = runif(n=1, min=0.5/2, max=0.5*1.5)
beta = runif(n=1, min=25/2, max=25*1.5)
gamma = runif(n=1, min=8/2, max=8*1.5)
sigmaW = runif(n=1, min=sqrt(10)/2, max=sqrt(10)*1.5)
sigmaV = runif(n=1, min=1/2, max=1*1.5)
# Initialize
initial.nlss <- NGMmodel(alpha, beta, gamma, sig2w=sigmaW^2, sig2v=sigmaV^2)
# Run
alltracesAve[,,irep] = CPFSAEM(initial.nlss=initial.nlss, y=obs,
n.particles=n.particles, n.em=n.em)$trace
cat("\n")
}
(time2run = proc.time() - ptm)
# graphics
truep = c(.5,25,8,sqrt(10),1)
names = c(expression(alpha), expression(beta), expression(gamma),
expression(sigma[w]), expression(sigma[v]))
dev.new(height=3)
par(mfcol=c(1,5), mar=c(3,1,1,1), mgp=c(1.6,.6,0), oma=c(0,1.5,1,0))
for (i in 1:5){
ts.plot(alltracesAve[,i,], ylab='', xlab='iteration')
mtext(names[i], line=.25, side=3, cex=1.1)
if (i == 1) mtext('trace', side=2, line=1.5, cex=.9)
abline(h=truep[i], col=2, lwd=1)
}
Example 12.8
library(plyr)
library(astsa)
data(jj)
y = jj
dev.new(height=5)
par(mar=c(2,2,0,0)+.7, mgp=c(1.6,.6,0))
plot(y, type='o', ylab='J&J QE/Share') # data shown in chapter 2
### setup - model and initial parameter values ###
# Model is written here as:
# x[t] = G x[t-1] + w[t]; w[t] ~ iid Np(0, W); p = 4
# y[t] = F' x[t] + v[t]; v[t] ~ iid Nq(0, V); q = 1
# x[0] ~ Np(a1, R1)
n = length(y)
F = c(1,1,0,0)
G = diag(0,4)
G[1,1] = 1.03 # this is φ
G[2,]=c(0,-1,-1,-1); G[3,]=c(0,1,0,0); G[4,]=c(0,0,1,0)
a1 = rbind(.7,0,0,0)
R1 = diag(.04,4)
V = .1
W11 = .1
W22 = .1
### FFBS scheme
ffbs = function(y,F,G,V,W11,W22,a1,R1){
n = length(y)
Ws = diag(c(W11,W22,1,1)) # 1s are added to ease coding- doesn't change the results
iW = diag(1/diag(Ws),4)
a = matrix(0,n,4)
R = array(0,c(n,4,4))
m = matrix(0,n,4)
C = array(0,c(n,4,4))
a[1,] = a1[,1]
R[1,,] = R1
f = t(F)%*%a[1,]
Q = t(F)%*%R[1,,]%*%F+V
A = R[1,,]%*%F/Q[1,1]
m[1,] = a[1,]+A%*%(y[1]-f)
C[1,,] = R[1,,]-A%*%t(A)*Q[1,1]
for (t in 2:n){
a[t,] = G%*%m[t-1,]
R[t,,] = G%*%C[t-1,,]%*%t(G) + Ws
f = t(F)%*%a[t,]
Q = t(F)%*%R[t,,]%*%F+V
A = R[t,,]%*%F/Q[1,1]
m[t,] = a[t,]+A%*%(y[t]-f)
C[t,,] = R[t,,]-A%*%t(A)*Q[1,1]
}
xb = matrix(0,n,4)
xb[n,] = m[n,] + t(chol(C[n,,]))%*%rnorm(4)
for (t in (n-1):1){
iC = solve(C[t,,])
CCC = solve(t(G)%*%iW%*%G+iC)
mmm = CCC%*%(t(G)%*%iW%*%xb[t+1,]+iC%*%m[t,])
xb[t,] = mmm + t(chol(CCC))%*%rnorm(4)
}
return(xb)
}
### MCMC scheme ###
#----------------
# Prior hyperparameters
# b0 = 0 # mean for beta = phi - 1
# B0 = Inf # var for beta (non-informative prior => use OLS for sampling beta)
n0 = 10 # use same for all variance components - these values get halved (1/2) below
s20v = .001 # for V
s20w =.05 # for Ws
set.seed(90210)
burnin = 1000
step = 20
M = 1000
niter = burnin+step*M
pars = matrix(0,niter,4)
xbs = array(0,c(niter,n,4))
pr <- progress_text() # displays progress (from plyr)
pr$init(niter)
for (iter in 1:niter){
xb = ffbs(y,F,G,V,W11,W22,a1,R1)
u = xb[,1] # trend components
yu = diff(u); xu = u[-n] # regress T[t]-T[t-1] on T[t-1]
regu = lm(yu~0+xu) # est of beta = phi - 1
phies = as.vector(coef(summary(regu)))[1:2] + c(1,0) # 1=estimate; 2=se
dft = df.residual(regu)
G[1,1] = phies[1] + rt(1,dft)*phies[2]
V = 1/rgamma(1, (n0+n)/2, (n0*s20v/2) + sum((y-xb[,1]-xb[,2])^2)/2)
W11 = 1/rgamma(1, (n0+n-1)/2, (n0*s20w/2) + sum((xb[-1,1]-phies[1]*xb[-n,1])^2)/2)
W22 = 1/rgamma(1, (n0+n-3)/2, (n0*s20w/2) + sum((xb[4:n,2] + xb[3:(n-1),2] +
xb[2:(n-2),2] +xb[1:(n-3),2])^2)/2)
xbs[iter,,] = xb
pars[iter,] = c(G[1,1], sqrt(V), sqrt(W11), sqrt(W22))
pr$step()
}
# Pix of Posteriors
ind = seq(burnin+1, niter, by=step)
names = c(expression(phi), expression(sigma[v]), expression(sigma[w~11]), expression(sigma[w~22]))
dev.new(height=6)
par(mfcol=c(3,4), mar=c(2,2,.25,0)+.75, mgp=c(1.6,.6,0), oma=c(0,0,1,0))
for (i in 1:4){
plot.ts(pars[ind,i], xlab="iterations", ylab="trace", main="")
mtext(names[i], side=3, line=.5, cex=1)
acf(pars[ind,i], main="", lag.max=25, xlim=c(1,25))
hist(pars[ind,i], main="", xlab="")
abline(v=mean(pars[ind,i]), lwd=2, col=3)
}
# Display smoothers of trend and season
mxb = cbind(apply(xbs[ind,,1],2,mean), apply(xbs[,,2],2,mean))
lxb = cbind(apply(xbs[ind,,1],2,quantile,0.005),apply(xbs[ind,,2],2,quantile,0.005))
uxb = cbind(apply(xbs[ind,,1],2,quantile,0.995),apply(xbs[ind,,2],2,quantile,0.995))
mxb=ts(mxb, start = tsp(jj)[1], freq=4)
lxb=ts(lxb, start = tsp(jj)[1], freq=4)
uxb=ts(uxb, start = tsp(jj)[1], freq=4)
names=c('Trend', 'Season')
dev.new(height=5)
par(mfrow=c(2,1), mar=c(2,2,0,0)+.7, mgp=c(1.6,.6,0))
for (i in 1:2){
L = min(lxb[,i])-.01; U = max(uxb[,i]) +.01
plot(mxb[,i], ylab=names[i], ylim=c(L,U), lwd=2)
xx=c(time(jj), rev(time(jj)))
yy=c(lxb[,i], rev(uxb[,i]))
polygon(xx, yy, border=NA, col=gray(.4, alpha = .3))
}
#########################################################################################
# This example and the next one are adapted from code by: Hedibert Freitas Lopes
#########################################################################################
Example 12.10
########################################################
## Some changes from the text - runs faster too ##
########################################################
# install.packages('plyr')
library(plyr) # used to view progress
## slight notation change in the code:
## state error st.dev. is 'sig' --- observation error st.dev. is 'tau'
# setup
mc.seed = 90210
burn.in = 1000
step.size = 3
n.mcmc = 2000*step.size + burn.in
g = function(arg){c(arg[1],arg[1]/(1+arg[1]^2),cos(1.2*arg[2]))}
fullxt = function(t,yt,xtm1,xt,xtp1,theta,sig,tau){
zm1 = g(c(xtm1,t-1))
z = g(c(xt,t))
full = dnorm(xt,sum(zm1*theta),tau,log=TRUE)+
dnorm(xtp1,sum(z*theta),tau,log=TRUE)+
dnorm(yt,xt^2/20,sig,log=TRUE)
return(full)
}
fullxn = function(t,yt,xtm1,xt,theta,sig,tau){
zm1 = g(c(xtm1,t-1))
z = g(c(xt,t))
full = dnorm(xt,sum(zm1*theta),tau,log=TRUE)+
dnorm(yt,xt^2/20,sig,log=TRUE)
return(full)
}
pr <- progress_text() # displays progress
pr$init(n.mcmc)
draw.x = function(y,x,x0,theta,sig,tau,v){
x1 = rnorm(1,x[1],v)
num = fullxt(1,y[1],x0,x1 ,x[2],theta,sig,tau)
den = fullxt(1,y[1],x0,x[1],x[2],theta,sig,tau)
if (log(runif(1))<(num-den)){x[1]=x1}
xn = rnorm(1,x[n],v)
num = fullxn(n,y[n],x[n-1],xn ,theta,sig,tau)
den = fullxn(n,y[n],x[n-1],x[n],theta,sig,tau)
if (log(runif(1))<(num-den)){x[n]=xn}
for (t in 2:(n-1)){
xt = rnorm(1,x[t],v)
num = fullxt(t,y[t],x[t-1],xt ,x[t+1],theta,sig,tau)
den = fullxt(t,y[t],x[t-1],x[t],x[t+1],theta,sig,tau)
if (log(runif(1))<(num-den)){x[t]=xt}
}
pr$step()
return(x)
}
fixedpar = function(y,X,b,A,v,lam){
n = length(y)
k = ncol(X)
par1 = v+n/2
var = solve(crossprod(X,X)+A)
mean = matrix(var%*%(crossprod(X,y)+crossprod(A,b)),k,1)
par2 = lam + sum((y-crossprod(t(X),mean))^2)
par2 = (par2 + crossprod(t(crossprod(mean-b,A)),mean-b))/2
sig2 = 1/rgamma(1,par1,par2)
var = var*sig2
mean = mean + crossprod(t(chol(var)),rnorm(k))
return(c(mean,sig2))
}
quant005 = function(x){quantile(x,0.005)}
quant995 = function(x){quantile(x,0.995)}
# Simulating the data
set.seed(mc.seed)
n = 100
sig2 = 1
tau2 = 10
sig = sqrt(sig2)
tau = sqrt(tau2)
theta = c(0.5,25,8)
ptr = c(theta,tau2,sig2)
y = rep(0,n)
x = rep(0,n)
x0 = 0
Z = g(c(x0,0))
x[1] = rnorm(1,sum(Z*theta),tau)
y[1] = rnorm(1,x[1]^2/20,sig)
for (t in 2:n){
Z = g(c(x[t-1],t-1))
x[t] = rnorm(1,sum(Z*theta),tau)
y[t] = rnorm(1,x[t]^2/20,sig)
}
xtr = x
# Prior hyperparameters
# note: IG prior originally was IG(a/2, ab/2), is now IG(a,b)
m0 = 0.0
C0 = 10.0
n0 = 6
sig20 = 5
theta0 = c(0.5,25,8)
V0 = diag(c(0.25,10,4))/tau2
nu0 = 6
tau20 = 50
iV0 = diag(1/diag(V0))
n0sig202 = sig20
# MCMC setup
v = .5 # tuning parameter
niter = n.mcmc
theta = ptr[1:3]
tau2 = ptr[4]
sig2 = ptr[5]
x = xtr
xs = NULL
ps = NULL
n0n2 = n0+n/2
ptm <- proc.time()
# Run it
for (iter in 1:niter){
x = draw.x(y,x,x0,theta,sig,tau,v)
sig2 = 1/rgamma(1,n0n2,n0sig202+sum((y-x^2/20)^2)/2)
sig = sqrt(sig2)
W = cbind(c(x0,x[1:(n-1)]),0:(n-1))
Z = t(apply(W,1,g))
par = fixedpar(x,Z,theta0,iV0,nu0,tau20)
theta = par[1:3]
tau2 = par[4]
tau = sqrt(tau2)
xs = rbind(xs,x)
ps = rbind(ps,c(par,sig2))
}
time2run = proc.time() - ptm
M0 = burn.in
M = n.mcmc-M0
L = step.size
index = seq((M0+1),niter,by=L)
nindex = length(index)
mx = apply(xs[index,],2,mean)
lx = apply(xs[index,],2,quant005)
ux = apply(xs[index,],2,quant995)
# display state estimation
dev.new(height=4)
par(mar=c(3,3,1.5,1), mgp=c(1.6,.6,0))
plot(xtr,type="p",pch=20,ylim=range(lx,ux+10,xtr),xlab="time",ylab="",main="State Estimation")
lines(mx,col=1,type="o")
legend("topleft",legend=c("True","Posterior Mean","99%CI"),col=c(1,1,gray(.5, alpha=.4)),
lty=c(0,1,1), pch=c(20,1,-1), lwd=c(1,1,5), cex=.8)
xx=c(time(xtr), rev(time(xtr)))
yy=c(lx, rev(ux))
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
# display parameter estimation
names = c(expression(alpha), expression(beta), expression(gamma),
expression(sigma[w]^2), expression(sigma[v]^2))
dev.new(height=5)
par(mfrow=c(3,5), mar=c(2.75,2.75,2,1), mgp=c(1.6,.6,0))
for (i in 1:5){
plot.ts(ps[index,i], xlab="iteration", ylab="")
title(names[i], cex=1.5)
abline(h=ptr[i],col=2,lwd=4)
}
for (i in 1:5)
acf(ps[index,i],main="", xlim=c(1.4,30))
for (i in 1:5){
hist(ps[index,i],prob=TRUE, xlab="",ylab="",main="")
lines(density(ps[index,i],adjust=2))
lp = apply(as.matrix(ps[index,i]),2,quant005)
up = apply(as.matrix(ps[index,i]),2,quant995)
abline(v=c(ptr[i],lp,up),col=c(2,4,4),lwd=c(4,1,1), lty=c(1,2,2))
}
time2run
Example 12.11
library(nltsa)
library(plyr)
npart = 1000 # Number of particles used in pmcmc
nmcmc = 1000 # Number of iterations in the mcmc samplers after burnin
burnin = 500 # Number of iterations to burn
step = 5 # Step size
mcmseed = 90210 # Beverly Hills seed
# Generate data
udat = NGMdata() # run at defaults [100 obs from the model given in the text]
y = udat$y # observations
X = udat$x # states
# Priors - IG(a,b) for variances and Normal(a,b) for regression coefs
sig2w00 = c(6, 50) # need a > 2 for variance ... or what kind of Bayesian are you?
sig2v00 = c(6, 5)
alpha00 = c(.5, sqrt(.25)) # note: means used for ALL parameters as initial value in pmcmc
beta00 = c(25, sqrt(10))
gamma00 = c(8, sqrt(4))
# Start it up
ptm <- proc.time() # measure how long this takes to run
u = pmcmcNGM(nmcmc, burnin, step, y, alpha00, beta00, gamma00, sig2w00, sig2v00, npart, mcmseed)
time2run = proc.time() - ptm
cat('acceptance rate =', u$accrate, '\n') # acceptance rate of the algorithm
# Some time in the future... we can make some nice pictures
parms = cbind(u$alpha, u$beta, u$gamma, u$sig2w, u$sig2v)
names = c(expression(alpha), expression(beta), expression(gamma),
expression(sigma[w]^2), expression(sigma[v]^2))
truep = c(.5, 25, 8, 10, 1)
dev.new(height=5.5)
par(mfcol=c(3,5), mar=c(3,3,1,1), mgp=c(1.6,.6,0), oma=c(0,0,1.5,0))
for (i in 1:5){
plot.ts(parms[,i], ylab='trace')
abline(h=truep[i], col=2, lwd=3)
mtext(names[i], 3, line=.5)
acf(parms[,i], 30, xlim=c(1.4,30)) # xlim lower bound is a cheesy way of hiding the 0 lag
hist(parms[,i], prob=TRUE, main="", xlab="")
lines(density(parms[,i], adjust=1.5))
lp = quantile(parms[,i], .005)
up = quantile(parms[,i], .995)
abline(v=c(truep[i],lp,up), col=c(2,4,4), lwd=c(3,1,1), lty=c(1,2,2))
}
dev.new(height=4)
par(mar=c(3,2,2,1), mgp=c(1.6,.6,0))
mX = apply(u$X, 2, mean)
lX = apply(u$X, 2, quantile,0.005)
uX = apply(u$X, 2, quantile,0.995)
plot.ts(mX, type='o', ylab="", main="State Estimation", ylim=c(-25,22))
xx=c(1:100, 100:1)
yy=c(lX, rev(uX))
polygon(xx, yy, border=NA, col=gray(.7, alpha = .4))
points(X, pch=20) # X is true state seq
legend("bottomleft",legend=c("True","Posterior Mean","99%CI"),col=c(4,1,gray(.6, alpha=.4)),
lty=c(0,1,1), pch=c(20,-1,-1), lwd=c(0,1,5), cex=.85)
###########################################################################
# The code for PMCMC here and PGAS in the next example is based on
# Matlab code from Fredrik Lindsten
###########################################################################
Example 12.12
library(nltsa)
library(plyr)
npart = 20 # Number of particles used in pmcmc
nmcmc = 2000 # Number of iterations in the mcmc samplers after burnin
burnin = 500 # Number of iterations to burn
mcmseed = 90210
# Generate data
udat = NGMdata() # run at defaults
y = udat$y # data
X = udat$x # states
# Priors - IG(a,b) for variance, Normal(mean,variance) for regression
sig2w00 = c(6, 50) # igamma parameters (a,b): mean=b/(a-1) used as initial value ...
sig2v00 = c(6, 5) # ... and need a > 2 for finite variance
alpha00 = c(.5, .25)
beta00 = c(25, 10)
gamma00 = c(8, 4)
# Run it and time it
ptm <- proc.time()
u = pgasNGM(nmcmc, burnin, y, alpha00, beta00, gamma00, sig2w00, sig2v00, npart, mcmseed)
(time2run = proc.time() - ptm)
# Pretty pictures
parms = cbind(u$theta[,1], u$theta[,2], u$theta[,3], u$sig2w, u$sig2v)
names = c(expression(alpha), expression(beta), expression(gamma),
expression(sigma[w]^2), expression(sigma[v]^2))
truep = c(.5, 25, 8, 10, 1)
dev.new(height=5.5)
par(mfcol=c(3,5), mar=c(3,3,1,1), mgp=c(1.6,.6,0), oma=c(0,0,1.5,0))
# Parameters ...
for (i in 1:5){
plot.ts(parms[,i], ylab='trace')
abline(h=truep[i], col=2, lwd=3)
mtext(names[i], 3, line=.25)
acf(parms[,i], 30, xlim=c(1.4,30))
hist(parms[,i], prob=TRUE, main="", xlab="")
lines(density(parms[,i], adjust=1.5))
lp = quantile(parms[,i], .005)
up = quantile(parms[,i], .995)
abline(v=c(truep[i],lp,up), col=c(2,4,4), lwd=c(3,1,1), lty=c(1,2,2))
}
# and States ...
dev.new(height=4)
par(mar=c(3,2,2,1), mgp=c(1.6,.6,0))
mX = apply(u$X, 2, mean)
lX = apply(u$X, 2, quantile,0.005)
uX = apply(u$X, 2, quantile,0.995)
plot.ts(mX, type='o', ylab="", main="State Estimation", ylim=c(min(lX,mX)-1,max(uX,mX)+8))
xx=c(1:100, 100:1)
yy=c(lX, rev(uX))
polygon(xx, yy, border=NA, col=gray(.7, alpha = .4))
points(X, pch=20) # X is true state seq
legend("topleft",legend=c("True","Posterior Mean","99%CI"),col=c(1,1,gray(.6, alpha=.4)),
lty=c(0,1,1), pch=c(20,1,-1), lwd=c(0,1,5), cex=.85)
