Nonlinear Time Series
Theory, Methods and Applications with R Examples
CRAN
Get R and related stuff
Home
Mr Jo
Randal's other job

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

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
rmspe = rep(0,n.ahead); x00 = ks\$xf[,,num]; P00 = ks\$Pf[,,num]
Q=t(cQ)%*%cQ;  R=t(cR)%*%(cR)
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))
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")
```

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="")
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="")
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="")
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)
```

© R.Douc, E.Moulines, D.Stoffer    Nonlinear Time Series: Theory, Methods and Applications with R Examples    Based on a design by: The All Art Directory