@@ -245,17 +245,17 @@ for (i in 1:5) {
245245 dens_conj_2 <- density(res_conj_2$sigma2, bw = "SJ", adj = 2)
246246 }
247247
248- plot(dens_improper, xlim = range(dens_improper$x, dens_semi_1$x, dens_semi_2$x) ,
249- ylim = range(dens_improper$y , dens_semi_1$y , dens_semi_2$y ),
250- main = "", xlab = name, ylab = "" )
248+ plot(dens_improper, main = "", xlab = name, ylab = "" ,
249+ xlim = range(dens_improper$x , dens_semi_1$x , dens_semi_2$x ),
250+ ylim = range(dens_improper$y, dens_semi_1$y, dens_semi_2$y) )
251251 if (i == 1) title("Semi-conjugate priors")
252252 lines(dens_semi_1, col = 2, lty = 2)
253253 lines(dens_semi_2, col = 3, lty = 3)
254254 legend("topright", c("Improper", "Tight", "Loose"), col = 1:3, lty = 1:3)
255255
256- plot(dens_improper, xlim = range(dens_improper$x, dens_conj_1$x, dens_conj_2$x) ,
257- ylim = range(dens_improper$y , dens_conj_1$y , dens_conj_2$y ),
258- main = "", xlab = name, ylab = "" )
256+ plot(dens_improper, main = "", xlab = name, ylab = "" ,
257+ xlim = range(dens_improper$x , dens_conj_1$x , dens_conj_2$x ),
258+ ylim = range(dens_improper$y, dens_conj_1$y, dens_conj_2$y) )
259259 if (i == 1) title("Conjugate priors")
260260 lines(dens_conj_1, col = 2, lty = 2)
261261 lines(dens_conj_2, col = 3, lty = 3)
@@ -361,7 +361,8 @@ To assess the probability of nonstationarity, we can simply count the
361361draws outside of the stationarity region.
362362
363363``` {r}
364- nonstationary <- matrix(NA, nrow(ar3draws), 3, dimnames = list(NULL, order = 1:3))
364+ nonstationary <- matrix(NA, nrow(ar3draws), 3,
365+ dimnames = list(NULL, order = 1:3))
365366nonstationary[, 1] <- abs(ar1draws) > 1
366367nonstationary[, 2] <- ar2draws[,1] + ar2draws[,2] > 1 |
367368 abs(ar2draws[,2]) > 1 |
@@ -572,8 +573,9 @@ for (p in 1:2) lines(density(mu[, as.character(p)], bw = "SJ", adj = 2),
572573 col = p + 1, lty = p + 1)
573574legend("topright", paste0("p = ", 0:2), col = 1:3, lty = 1:3)
574575
575- plot(density(sigma2[, "0"], bw = "SJ", adj = 2), xlim = range(sigma2), ylab = "",
576- main = "Posterior of the marginal variance", xlab = expression(sigma^2))
576+ plot(density(sigma2[, "0"], bw = "SJ", adj = 2), xlim = range(sigma2),
577+ xlab = expression(sigma^2), ylab = "",
578+ main = "Posterior of the marginal variance")
577579for (p in 1:2) lines(density(sigma2[, as.character(p)], bw = "SJ", adj = 2),
578580 col = p + 1, lty = p + 1)
579581legend("topright", paste0("p = ", 0:2), col = 1:3, lty = 1:3)
@@ -621,8 +623,9 @@ for (p in 1:2) lines(density(mu[[as.character(p)]], bw = "SJ", adj = 2),
621623 col = p + 1, lty = p + 1)
622624legend("topright", paste0("p = ", 0:2), col = 1:3, lty = 1:3)
623625
624- plot(density(sigma2[["0"]], bw = "SJ", adj = 2), xlim = range(sigma2), ylab = "",
625- main = "Posterior of the marginal variance", xlab = expression(sigma^2))
626+ plot(density(sigma2[["0"]], bw = "SJ", adj = 2), xlim = range(sigma2),
627+ xlab = expression(sigma^2), ylab = "",
628+ main = "Posterior of the marginal variance", )
626629for (p in 1:2) lines(density(sigma2[[as.character(p)]], bw = "SJ", adj = 2),
627630 col = p + 1, lty = p + 1)
628631legend("topright", paste0("p = ", 0:2), col = 1:3, lty = 1:3)
@@ -700,7 +703,8 @@ for (m in seq_len(ndraws + nburn)) {
700703 if (-1 < phiprop & phiprop < 1) {
701704 logR <- dbetarescaled(phiprop, aphi, bphi, log = TRUE) -
702705 dbetarescaled(phi, aphi, bphi, log = TRUE) +
703- dnorm(y0, zeta / (1 - phiprop), sqrt(sigma2 / (1 - phiprop^2)), log = TRUE) -
706+ dnorm(y0, zeta / (1 - phiprop), sqrt(sigma2 / (1 - phiprop^2)),
707+ log = TRUE) -
704708 dnorm(y0, zeta / (1 - phi), sqrt(sigma2 / (1 - phi^2)), log = TRUE)
705709 if (log(runif(1)) < logR) {
706710 phi <- phiprop
@@ -769,12 +773,14 @@ for (m in seq_len(ndraws + nburn)) {
769773 # Step (d): Draw the persistence
770774 tmp <- y0^2 + sum(y[-length(y)]^2)
771775 propvar <- 1 / (1 / priorvar + tmp / sigma2)
772- propmean <- propvar * (priormean / priorvar + (y0 * (y[1] - zeta) + sum(y[-length(y)] * (y[-1] - zeta))) / sigma2)
776+ propmean <- propvar * (priormean / priorvar + (y0 * (y[1] - zeta) +
777+ sum(y[-length(y)] * (y[-1] - zeta))) / sigma2)
773778 phiprop <- rnorm(1, propmean, sqrt(propvar))
774779 if (-1 < phiprop & phiprop < 1) {
775780 logR <- dbetarescaled(phiprop, aphi, bphi, log = TRUE) -
776781 dbetarescaled(phi, aphi, bphi, log = TRUE) +
777- dnorm(y0, zeta / (1 - phiprop), sqrt(sigma2 / (1 - phiprop^2)), log = TRUE) -
782+ dnorm(y0, zeta / (1 - phiprop), sqrt(sigma2 / (1 - phiprop^2)),
783+ log = TRUE) -
778784 dnorm(y0, zeta / (1 - phi), sqrt(sigma2 / (1 - phi^2)), log = TRUE) +
779785 dnorm(phi, priormean, sqrt(priorvar), log = TRUE) -
780786 dnorm(phiprop, priormean, sqrt(priorvar), log = TRUE)
@@ -872,12 +878,14 @@ for (m in seq_len(ndraws + nburn)) {
872878 # Step (d): Draw the persistence
873879 tmp <- y0^2 + sum(y[-length(y)]^2)
874880 propvar <- 1 / (1 / priorvar + tmp / sigma2)
875- propmean <- propvar * (priormean / priorvar + (y0 * (y[1] - zeta) + sum(y[-length(y)] * (y[-1] - zeta))) / sigma2)
881+ propmean <- propvar * (priormean / priorvar +
882+ (y0 * (y[1] - zeta) + sum(y[-length(y)] * (y[-1] - zeta))) / sigma2)
876883 phiprop <- rnorm(1, propmean, sqrt(propvar))
877884 if (-1 < phiprop & phiprop < 1) {
878885 logR <- dbetarescaled(phiprop, aphi, bphi, log = TRUE) -
879886 dbetarescaled(phi, aphi, bphi, log = TRUE) +
880- dnorm(y0, zeta / (1 - phiprop), sqrt(sigma2 / (1 - phiprop^2)), log = TRUE) -
887+ dnorm(y0, zeta / (1 - phiprop), sqrt(sigma2 / (1 - phiprop^2)),
888+ log = TRUE) -
881889 dnorm(y0, zeta / (1 - phi), sqrt(sigma2 / (1 - phi^2)), log = TRUE) +
882890 dnorm(phi, priormean, sqrt(priorvar), log = TRUE) -
883891 dnorm(phiprop, priormean, sqrt(priorvar), log = TRUE)
@@ -1039,13 +1047,13 @@ matrices with values $N_{i,hk}$.
10391047getTransitions <- function(x, classes) {
10401048 transitions <- matrix(0, length(classes), length(classes))
10411049 for (i in seq_len(length(x) - 1)) {
1042- transitions[x[i], x[i+ 1]] <- transitions[x[i], x[i+ 1]] + 1
1050+ transitions[x[i], x[i + 1]] <- transitions[x[i], x[i + 1]] + 1
10431051 }
10441052 dimnames(transitions) <- list(from = classes, to = classes)
10451053 transitions
10461054}
10471055income_transitions <- lapply(seq_len(nrow(income)),
1048- function(i) getTransitions(income[i,], classes = 0:5))
1056+ function(i) getTransitions(income[i,], classes = 0:5))
10491057```
10501058
10511059Based on the transition matrices, the total transitions between the
@@ -1061,7 +1069,8 @@ knitr::kable(income_trans_male)
10611069We obtain the posterior mean estimates based on a uniform prior.
10621070
10631071``` {r}
1064- income_trans_female <- (1 + income_trans_female) / rowSums(1 + income_trans_female)
1072+ income_trans_female <- (1 + income_trans_female) /
1073+ rowSums(1 + income_trans_female)
10651074income_trans_male <- (1 + income_trans_male) / rowSums(1 + income_trans_male)
10661075knitr::kable(income_trans_female, digits = 3)
10671076knitr::kable(income_trans_male, digits = 3)
@@ -1073,9 +1082,8 @@ if (pdfplots) {
10731082 par(mar = c(2.6, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 2)
10741083}
10751084par(mfrow = c(1, 2))
1076- library("corrplot")
1077- corrplot(income_trans_female, method = "square",
1078- is.corr = FALSE, col = 1, cl.pos = "n")
1079- corrplot(income_trans_male, method = "square",
1080- is.corr = FALSE, col = 1, cl.pos = "n")
1085+ corrplot::corrplot(income_trans_female, method = "square", is.corr = FALSE,
1086+ col = 1, cl.pos = "n")
1087+ corrplot::corrplot(income_trans_male, method = "square", is.corr = FALSE,
1088+ col = 1, cl.pos = "n")
10811089```
0 commit comments