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

R Code

Click [+] to expand or collapse section. Click [-] to collapse entire page. Or [+ Expand the Entire Page +]

Before you start, open an R session, download nltsa ONCE as follows:

   ######################## do this only once ############################
   ##                                                                   ##       
download.file("http://www.stat.pitt.edu/stoffer/nltsa/nltsa.rda", "nltsa.rda")
   ##                                                                   ## 
   ################# copy and paste the above line into R ################ 

Then, like a package, when needed (but 'load' instead of 'require')

 load('nltsa.rda')

If you save the workspace when you're done, you don't have to keep loading the file. Otherwise, load nltsa.rda each time you use it (you only download it once). There is sort of a NLTSA MANUAL, but it's not very useful without the text. When you see require('...') or library('...'), it is assumed you have installed the corresponding package.





[+] Chapter 1

    Example 1.36

    require(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

    require(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

    require(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

    require(astsa)
    require(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

    require(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

    require(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

    require(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

    require(astsa)
    data(flu)
    
    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

    require(astsa)
    data(EXP6)
    
    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

    require(vcd)                # load package - install it if you don't have it
    require(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

    require(tsDyn)       # load package - install it if you don't have it
    require(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

    require(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 --## 
    require(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 --## 
    require(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

    load('nltsa.rda')
    
    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

    require(TSA)    # load package - install it if you don't have it
    Keenan.test(log10(lynx))
    Tsay.test(log10(lynx))
    

[+] Chapter 4

    Figure 4.1

    require(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

    require(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

    load('nltsa.rda')
    
    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

    load('nltsa.rda')
    
    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
    require(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

    require(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

    require(astsa)
    require(bayesGARCH)   # load package - install it if you don't have it
    require(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

    require(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
    # require(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

    load('nltsa.rda')
    
    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

    require(depmixS4)
    load('nltsa.rda')
    
    ## 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

    require(depmixS4)
    load('nltsa.rda')
    
    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

    require(astsa)
    data(flu)
    require(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)
    load('nltsa.rda') 
    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)
    load('nltsa.rda')
    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)
    load('nltsa.rda')
    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)
    load('nltsa.rda')
    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)
    load('nltsa.rda')
    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)
    load('nltsa.rda')
    
    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(ggplot2)
    library(plyr)
    library(reshape2)
    library(testthat)
    load('nltsa.rda')
    
    #######################
    ##### 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)))
    
    # 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,,])
    expect_true(all(smoothed.ffbsi$weights == 1))
    expect_true(all(smoothed.ffbsilin$weights == 1))
    
    # 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
    
    # 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.df <- rbind(traj.poor.df,
                     traj.ffbsm.df,
                     traj.ffbsi.df,
                     traj.ffbsilin.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"))
    
    dev.new(height=5)
    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"))
    
    #======= Study the variance of the posterior-mean estimate ===========
    # Compute nrep posterior mean estimates: a bit long to run
    estimates.poor <- array(data=NA_real_, dim=c(nrep,T),dimnames=c("Repetition","Time"))
    estimates.ffbsm <- array(data=NA_real_, dim=c(nrep,T),dimnames=c("Repetition","Time"))
    estimates.ffbsi <- array(data=NA_real_, dim=c(nrep,T),dimnames=c("Repetition","Time"))
    estimates.ffbsilin <- array(data=NA_real_, dim=c(nrep,T),dimnames=c("Repetition","Time"))
    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)))
      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,,])
        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.df <- rbind(est.poor.df,
                    est.ffbsm.df,
                    est.ffbsi.df,
                    est.ffbsilin.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"))
    
    dev.new(height=5)
    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

    load('nltsa.rda')    
    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

    load('nltsa.rda')
    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

    require(plyr)
    require(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

    load('nltsa.rda')
    require(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    = 25        #  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)
    		  
    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

    load('nltsa.rda')
    require(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)
    



  © 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