This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* ----------------------------------------------------------------------------- | |
Estimate latent variable hierarchical Bayes RFM model | |
stan ver. 2.14 | |
CREATED: by ill-identified, 18/APR/2016 | |
REVISED: 25/FEB/2017, REVISED | |
------------------------------------------------------------------------------- */ | |
data { |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(dplyr) | |
library(tidyr) | |
library(ggplot2) | |
library(rstan) | |
library(loo) | |
library(ggmcmc) | |
# read datasets | |
df <- read.csv("rfm.csv", stringsAsFactors = F) | |
colnames(df)[1] <- "ID" |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(ggmcmc) # ver. 1.1 | |
library(data.table) | |
ggs2 <- function(result, pattern=NULL, invert=F, ...){ | |
if(class(result)[1] =="stanfit" ){ | |
nThin <- result@sim$thin | |
nChains <- result@sim$chains | |
nBurnin <- result@sim$warmup2[1] | |
nIterations <- result@sim$n_save[1] - nBurnin | |
result <- as.data.frame(result) %>% | |
mutate(Chain=rep(1:nChains, nIterations), Iteration=rep(1:nIterations, nChains)) %>% |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
glmnet.mock <- function(formula=y~., family="gaussian", data, ...){ | |
fam.link <- strsplit(family, "\\(|\\)")[[1]] | |
family <- fam.link[1] | |
if(length(fam.link) >=2) link <- fam.link[2] | |
else link <- "identity" | |
return(glmnet(x=model.matrix(formula, data), | |
y=get(link)(model.response(model.frame(formula, data))), | |
family = family, | |
...) | |
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
require(mlogit) # 0.2-4 | |
as.mldata <- function(data){ | |
# convert HC dataset | |
# The alternatives are | |
# 1. Gas central heat with cooling (gcc) | |
# 2. Electric central resistence heat with cooling (ecc) | |
# 3. Electric room resistence heat with cooling (erc) | |
# 4. Electric heat pump, which provides cooling also (hpc) | |
# 5. Gas central heat without cooling (gc) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
source("common.R", encoding = "utf-8") | |
# ---- dlm ------ | |
# model 3 | |
res <- data.frame() | |
for( a in seq(from=.1, to=.95, by=.05) ){ | |
RP <- calc_RP(RP = df$RP095, p = df$AP, a = a) | |
z <- calc_Z(RP = RP, p = df$AP) | |
z1 <- z$z1 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
source("common.R", encoding = "utf-8") | |
# ----- KFAS ------ | |
df$RP <- calc_RP(df$RP095, df$AP, .95) | |
z <- calc_Z(RP = df$RP, p = df$AP) | |
df <- mutate(df, z1=z$z1, z2=z$z2) | |
df <- mutate(df, z1E=z1*end, z2E=z2*end) | |
# specify model | |
model3KFAS <- SSModel(logPI ~ SSMtrend(1, Q=NA) + |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
source("common.R", encoding = "utf-8") | |
df$RP <- calc_RP(df$RP095, df$AP, .95) | |
z <- calc_Z(RP = df$RP, p = df$AP) | |
df <- mutate(df, z1=z$z1, z2=z$z2) | |
df <- mutate(df, z1E=z1*end, z2E=z2*end) | |
ss <- AddLocalLevel(list(), y = df$logPI) # c | |
ss <- AddAr(ss, lags=2, y = df$logPI) # AR(2) | |
# time-varying regression |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# ------ common part ---- | |
require(ggplot2) | |
require(dplyr) | |
require(tidyr) | |
require(dlm) # 1.1-4 | |
require(KFAS) # 1.2.9 | |
require(bsts) # 0.7.1 | |
# calculate the reference price |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
require(bsts) # 0.7.1 | |
data(iclaims) # bring the initial.claims data into scope | |
# --- model 1 ---- | |
ss <- AddLocalLinearTrend(list(), initial.claims$iclaimsNSA) | |
ss <- AddSeasonal(ss, initial.claims$iclaimsNSA, nseasons = 52) | |
model1 <- bsts(initial.claims$iclaimsNSA, | |
state.specification = ss, | |
niter = 1000) |