Include libraries
set.seed(111)
library(readr)
## Warning: パッケージ 'readr' はバージョン 4.3.3 の R の下で造られました
library(ggplot2)
## Warning: パッケージ 'ggplot2' はバージョン 4.3.3 の R の下で造られました
library(rstan)
## 要求されたパッケージ StanHeaders をロード中です
## Warning: パッケージ 'StanHeaders' はバージョン 4.3.3 の R の下で造られました
##
## rstan version 2.32.6 (Stan version 2.32.2)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
library(rstudioapi)
## Warning: パッケージ 'rstudioapi' はバージョン 4.3.3 の R の下で造られました
Read data
m <- read_csv("leadingIndex.csv", show_col_types=FALSE)
# Composite Leading Index Data
# Date,LeadingIndex
# 1985-01-01,85.5
# 1985-02-01,85.7
# ...............
# Peaks and Troughs for business cycles
recess.df = read.table(textConnection(
"Peak, Trough
1985-06-01, 1986-11-01
1991-02-01, 1993-10-01
1997-05-01, 1999-01-01
2000-11-01, 2002-01-01
2008-02-01, 2009-03-01
2012-03-01, 2012-11-01
2018-10-01, 2020-05-01 "), sep=',',
colClasses=c('Date', 'Date'), header=TRUE)
Time series plot of Composite Index (shaded area: recession)
month <- as.vector(t(m[,1]))
ci <- as.vector(t(m[,2]))
#
ci.df <-data.frame(dates = as.Date(month), value = ci)
T <- length(ci)
#
g1 <- ggplot(ci.df) + geom_line(aes(x=dates, y=value)) +
xlab("Time") + ylab("Composite Index")
g1 = g1 + geom_rect(data=recess.df,aes(xmin=Peak, xmax=Trough, ymin=-Inf, ymax=+Inf), fill = "red", alpha=0.2)
g1
Time series plot of Change in Composite Index (%)
y <- ( log(ci[2:T]) - log(ci[1:T-1]) ) * 100
month <- month[2:T]
T <- T-1
cirate.df <-data.frame(dates = as.Date(month), value = y)
#
g2 <- ggplot(cirate.df) + geom_line(aes(x=dates, y=value)) +
xlab("Time") + ylab("Composite Index Change")
g2 = g2 + geom_rect(data=recess.df,aes(xmin=Peak, xmax=Trough, ymin=-Inf, ymax=+Inf), fill = "red", alpha=0.2)
g2
Markov Switching Model using rstan
mydata <- list( T = T, y = y);
myinit <- function() {
list( mu1 = -1, mu2 = 1, sigma2 = 1, p11 = 0.90, p22 = 0.95)
}
fit <- stan(file = 'markov.stan', data = mydata, init = myinit,
chains = 4, warmup = 1000, iter = 11000, thin=10, cores = parallel::detectCores())
Estimation results
print(fit, pars = c("mu1","mu2","sigma2","p11","p22","loglik"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=11000; warmup=1000; thin=10;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu1 -0.37 0.00 0.13 -0.67 -0.44 -0.34 -0.27 -0.17 3987 1
## mu2 0.92 0.00 0.25 0.48 0.73 0.95 1.11 1.35 4173 1
## sigma2 0.59 0.00 0.06 0.48 0.55 0.59 0.63 0.72 4094 1
## p11 0.93 0.00 0.02 0.88 0.92 0.93 0.94 0.96 4045 1
## p22 0.92 0.00 0.02 0.87 0.90 0.92 0.93 0.96 4013 1
## loglik -745.44 0.03 2.04 -749.82 -746.77 -745.30 -743.86 -742.14 3914 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 6 14:21:09 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
param <-c("mu1","mu2","sigma2","p11","p22")
stan_dens(fit, pars = param, nrow = 3, ncol = 2, separate_chains = TRUE)
stan_trace(fit, pars = param, nrow = 3, ncol = 2)
stan_ac(fit, pars = param, nrow = 3, ncol = 2, separate_chains = TRUE)
AIC, BIC. # param = 5 (mu1, mu2, sigma2, p11, p22)
mcmcout = rstan::extract(fit)
cat("AIC=",-2*max(mcmcout$loglik)+2*5,"\n")
## AIC= 1492.467
cat("BIC=",-2*max(mcmcout$loglik)+log(T)*5)
## BIC= 1513.198
