Time series fit and posterior predictive model checking for unemployment series. See Chapter 11 in Regression and Other Stories.


Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("ggplot2")
library("bayesplot")
theme_set(bayesplot::theme_default(base_family = "sans"))

Load data

unemp <- read.table(root("Unemployment/data","unemp.txt"), header=TRUE)
head(unemp)
  year   y
1 1947 3.9
2 1948 3.8
3 1949 5.9
4 1950 5.3
5 1951 3.3
6 1952 3.0

Plot the unemployment rate

par(mar=c(3,3,1,.1), mgp=c(1.7,.5,0), tck=-.01)
plot(unemp$year, unemp$y, type="l", ylab="Unemployment rate", xlab="Year", yaxs="i",
  ylim=c(0, max(unemp$y)*1.05), xaxt="n", yaxt="n", bty="l")
axis(1, seq(1950,2010,10), rep("",7))
axis(1, seq(1950,2010,20))
axis(2, seq(0,10), rep("",11))
axis(2, c(0,5,10), paste (c(0,5,10), "%", sep=""))

Fit a 1st-order autogregression

n <- nrow(unemp)
unemp$y_lag <- c(NA, unemp$y[1:(n-1)])
fit_lag <- stan_glm(y ~ y_lag, data=unemp, refresh=0)
print(fit_lag, digits=2)
stan_glm
 family:       gaussian [identity]
 formula:      y ~ y_lag
 observations: 69
 predictors:   2
------
            Median MAD_SD
(Intercept) 1.35   0.45  
y_lag       0.77   0.08  

Auxiliary parameter(s):
      Median MAD_SD
sigma 1.03   0.09  

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Simulate replicated datasets using posterior predict

y_rep <- posterior_predict(fit_lag)
y_rep <- cbind(unemp$y[1], y_rep)
n_sims <- nrow(y_rep)

Simulate replicated datasets "manually"

sims <- as.matrix(fit_lag)
n_sims <- nrow(sims)
y_rep <- array(NA, c(n_sims, n))
for (s in 1:n_sims){
  y_rep[s,1] <- unemp$y[1]
  for (t in 2:n){
    y_rep[s,t] <- sims[s,"(Intercept)"] + sims[s,"y_lag"] * y_rep[s,t-1] + rnorm(1, 0, sims[s,"sigma"])
  }
}

Plot the simulated unemployment rate series

par(mar=c(1,1,3,.1), mgp=c(2,.5,0), tck=-.01)
par(mfrow=c(3,5))
for (s in sort(sample(n_sims, 15))){
  plot (unemp$year, y_rep[s,], type="l", ylab="", xlab="", yaxs="i",
  ylim=c(0, max(unemp$y)*1.05), xaxt="n", yaxt="n", bty="l", main=paste("Simulation", s))
  axis(1, seq(1950,2010,10), rep("",7))
  axis(2, seq(0,10), rep("",11))
}

Numerical posterior predictive check

test <- function (y){
  n <- length(y)
  y_lag <- c(NA, y[1:(n-1)])
  y_lag_2 <- c(NA, NA, y[1:(n-2)])
  return(sum(sign(y-y_lag) != sign(y_lag-y_lag_2), na.rm=TRUE))
}
test_y <- test(unemp$y)
test_rep <- apply(y_rep, 1, test)
print(mean(test_rep > test_y))
[1] 0.993
print(quantile(test_rep, c(.1,.5,.9)))
10% 50% 90% 
 31  36  42 

Plot test statistic for data and histogram of test statistics for replications

ppc_stat(y=unemp$y, yrep=y_rep, stat=test, binwidth = 1) + scale_y_continuous(breaks=NULL)

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogVW5lbXBsb3ltZW50IgphdXRob3I6ICJBbmRyZXcgR2VsbWFuLCBKZW5uaWZlciBIaWxsLCBBa2kgVmVodGFyaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICB0b2M6IHRydWUKICAgIHRvY19kZXB0aDogMgogICAgdG9jX2Zsb2F0OiB0cnVlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQpUaW1lIHNlcmllcyBmaXQgYW5kIHBvc3RlcmlvciBwcmVkaWN0aXZlIG1vZGVsIGNoZWNraW5nIGZvcgp1bmVtcGxveW1lbnQgc2VyaWVzLiBTZWUgQ2hhcHRlciAxMSBpbiBSZWdyZXNzaW9uIGFuZCBPdGhlcgpTdG9yaWVzLgoKLS0tLS0tLS0tLS0tLQoKCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGNvbW1lbnQ9TkEpCiMgc3dpdGNoIHRoaXMgdG8gVFJVRSB0byBzYXZlIGZpZ3VyZXMgaW4gc2VwYXJhdGUgZmlsZXMKc2F2ZWZpZ3MgPC0gRkFMU0UKYGBgCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyIH0KbGlicmFyeSgicnByb2pyb290IikKcm9vdDwtaGFzX2ZpbGUoIi5ST1MtRXhhbXBsZXMtcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpsaWJyYXJ5KCJyc3RhbmFybSIpCmxpYnJhcnkoImdncGxvdDIiKQpsaWJyYXJ5KCJiYXllc3Bsb3QiKQp0aGVtZV9zZXQoYmF5ZXNwbG90Ojp0aGVtZV9kZWZhdWx0KGJhc2VfZmFtaWx5ID0gInNhbnMiKSkKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CiMgZ3JheXNjYWxlIGZpZ3VyZXMgZm9yIHRoZSBib29rCmlmIChzYXZlZmlncykgY29sb3Jfc2NoZW1lX3NldChzY2hlbWUgPSAiZ3JheSIpCmBgYAoKIyMjIyBMb2FkIGRhdGEKCmBgYHtyIH0KdW5lbXAgPC0gcmVhZC50YWJsZShyb290KCJVbmVtcGxveW1lbnQvZGF0YSIsInVuZW1wLnR4dCIpLCBoZWFkZXI9VFJVRSkKaGVhZCh1bmVtcCkKYGBgCgojIyMjIFBsb3QgdGhlIHVuZW1wbG95bWVudCByYXRlCgpgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJVbmVtcGxveW1lbnQvZmlncyIsInVuZW1wbG95bWVudDEucGRmIiksIGhlaWdodD0zLCB3aWR0aD00LjUpCmBgYApgYGB7ciB9CnBhcihtYXI9YygzLDMsMSwuMSksIG1ncD1jKDEuNywuNSwwKSwgdGNrPS0uMDEpCnBsb3QodW5lbXAkeWVhciwgdW5lbXAkeSwgdHlwZT0ibCIsIHlsYWI9IlVuZW1wbG95bWVudCByYXRlIiwgeGxhYj0iWWVhciIsIHlheHM9ImkiLAogIHlsaW09YygwLCBtYXgodW5lbXAkeSkqMS4wNSksIHhheHQ9Im4iLCB5YXh0PSJuIiwgYnR5PSJsIikKYXhpcygxLCBzZXEoMTk1MCwyMDEwLDEwKSwgcmVwKCIiLDcpKQpheGlzKDEsIHNlcSgxOTUwLDIwMTAsMjApKQpheGlzKDIsIHNlcSgwLDEwKSwgcmVwKCIiLDExKSkKYXhpcygyLCBjKDAsNSwxMCksIHBhc3RlIChjKDAsNSwxMCksICIlIiwgc2VwPSIiKSkKYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZGV2Lm9mZigpCmBgYAoKIyMjIyBGaXQgYSAxc3Qtb3JkZXIgYXV0b2dyZWdyZXNzaW9uCgpgYGB7ciB9Cm4gPC0gbnJvdyh1bmVtcCkKdW5lbXAkeV9sYWcgPC0gYyhOQSwgdW5lbXAkeVsxOihuLTEpXSkKZml0X2xhZyA8LSBzdGFuX2dsbSh5IH4geV9sYWcsIGRhdGE9dW5lbXAsIHJlZnJlc2g9MCkKcHJpbnQoZml0X2xhZywgZGlnaXRzPTIpCmBgYAoKIyMjIyBTaW11bGF0ZSByZXBsaWNhdGVkIGRhdGFzZXRzIHVzaW5nIHBvc3RlcmlvciBwcmVkaWN0CgpgYGB7ciB9CnlfcmVwIDwtIHBvc3Rlcmlvcl9wcmVkaWN0KGZpdF9sYWcpCnlfcmVwIDwtIGNiaW5kKHVuZW1wJHlbMV0sIHlfcmVwKQpuX3NpbXMgPC0gbnJvdyh5X3JlcCkKYGBgCgojIyMjIFNpbXVsYXRlIHJlcGxpY2F0ZWQgZGF0YXNldHMgIm1hbnVhbGx5IgoKYGBge3IgfQpzaW1zIDwtIGFzLm1hdHJpeChmaXRfbGFnKQpuX3NpbXMgPC0gbnJvdyhzaW1zKQp5X3JlcCA8LSBhcnJheShOQSwgYyhuX3NpbXMsIG4pKQpmb3IgKHMgaW4gMTpuX3NpbXMpewogIHlfcmVwW3MsMV0gPC0gdW5lbXAkeVsxXQogIGZvciAodCBpbiAyOm4pewogICAgeV9yZXBbcyx0XSA8LSBzaW1zW3MsIihJbnRlcmNlcHQpIl0gKyBzaW1zW3MsInlfbGFnIl0gKiB5X3JlcFtzLHQtMV0gKyBybm9ybSgxLCAwLCBzaW1zW3MsInNpZ21hIl0pCiAgfQp9CmBgYAoKIyMjIyBQbG90IHRoZSBzaW11bGF0ZWQgdW5lbXBsb3ltZW50IHJhdGUgc2VyaWVzCgpgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQppZiAoc2F2ZWZpZ3MpIHBkZihyb290KCJVbmVtcGxveW1lbnQvZmlncyIsInVuZW1wbG95bWVudDIucGRmIiksIGhlaWdodD00LjUsIHdpZHRoPTcuNSkKYGBgCmBgYHtyIH0KcGFyKG1hcj1jKDEsMSwzLC4xKSwgbWdwPWMoMiwuNSwwKSwgdGNrPS0uMDEpCnBhcihtZnJvdz1jKDMsNSkpCmZvciAocyBpbiBzb3J0KHNhbXBsZShuX3NpbXMsIDE1KSkpewogIHBsb3QgKHVuZW1wJHllYXIsIHlfcmVwW3MsXSwgdHlwZT0ibCIsIHlsYWI9IiIsIHhsYWI9IiIsIHlheHM9ImkiLAogIHlsaW09YygwLCBtYXgodW5lbXAkeSkqMS4wNSksIHhheHQ9Im4iLCB5YXh0PSJuIiwgYnR5PSJsIiwgbWFpbj1wYXN0ZSgiU2ltdWxhdGlvbiIsIHMpKQogIGF4aXMoMSwgc2VxKDE5NTAsMjAxMCwxMCksIHJlcCgiIiw3KSkKICBheGlzKDIsIHNlcSgwLDEwKSwgcmVwKCIiLDExKSkKfQpgYGAKYGBge3IgZXZhbD1GQUxTRSwgaW5jbHVkZT1GQUxTRX0KaWYgKHNhdmVmaWdzKSBkZXYub2ZmKCkKYGBgCgojIyMjIE51bWVyaWNhbCBwb3N0ZXJpb3IgcHJlZGljdGl2ZSBjaGVjawoKYGBge3IgfQp0ZXN0IDwtIGZ1bmN0aW9uICh5KXsKICBuIDwtIGxlbmd0aCh5KQogIHlfbGFnIDwtIGMoTkEsIHlbMToobi0xKV0pCiAgeV9sYWdfMiA8LSBjKE5BLCBOQSwgeVsxOihuLTIpXSkKICByZXR1cm4oc3VtKHNpZ24oeS15X2xhZykgIT0gc2lnbih5X2xhZy15X2xhZ18yKSwgbmEucm09VFJVRSkpCn0KdGVzdF95IDwtIHRlc3QodW5lbXAkeSkKdGVzdF9yZXAgPC0gYXBwbHkoeV9yZXAsIDEsIHRlc3QpCnByaW50KG1lYW4odGVzdF9yZXAgPiB0ZXN0X3kpKQpwcmludChxdWFudGlsZSh0ZXN0X3JlcCwgYyguMSwuNSwuOSkpKQpgYGAKCiMjIyMgUGxvdCB0ZXN0IHN0YXRpc3RpYyBmb3IgZGF0YSBhbmQgaGlzdG9ncmFtIG9mIHRlc3Qgc3RhdGlzdGljcyBmb3IgcmVwbGljYXRpb25zCgpgYGB7ciB9CnBwY19zdGF0KHk9dW5lbXAkeSwgeXJlcD15X3JlcCwgc3RhdD10ZXN0LCBiaW53aWR0aCA9IDEpICsgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcz1OVUxMKQpgYGAKCg==