generate.Rmd

---
title: "Generate Data and Test"
author: "Amos B. Elberg"
date: "2/7/2019"
output: html_document
params:
  ####### Parameters on data generation
  P: !r 50                  # Number of periods
  use_sampled: F                  # Whether to use the exact prices generated or sample
  N: 700                          # Number of sampled prices if use_sampled = T
  S: 5                           # Number of price series
  F: 2                            # Number of features for regression
  ###### Parameters controlling the model 
  periods_to_predict: 30
  period_scale: 52                # Prior on the scale of price movements
  # a reasonable starting place is the number of periods over which 
  # its reasonably expected that a price may double; e.g., 52 for weekly prices
  p: 3                            
  q: 3
  ar: 3                           # Number of periods for trend model
  s: !r c(4, 12)                  # Seasonality periods
  corr_prior: 4
  ###### Parameters controlling hierarchical shrinkage 
  nu: 3
  m0: .05                         # Prior on proportion of coefficients that are relevant
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(magrittr)
library(ggplot2)
library(rethinking)
library(bayesplot)
library(tidybayes)
library(rstan)
library(ggplot2)
library(forcats)
library(purrr)
library(rethinking)
library(MASS)
library(tidyr)
library(ccgarch)
library(abind)
library(marima)
```

# Generate Fake Data 

```{r seed}
set.seed(1974)
periods_to_generate <- params$P + params$periods_to_predict
period_vline <- data.frame(xintercept=params$P)
```

## Trends

```{r trends}
theta_trend <- abs(rcauchy(params$S, 0, .01)) / params$period_scale
omega_trend <- rlkjcorr(1, params$S, eta=params$corr_prior)
alpha_trend <- rnorm(params$S, 0, .01) / params$period_scale
#beta_trend <- pmin(abs(matrix(rnorm(params$S * params$ar, sd=0.1),  ncol=params$S)), .3) 
L_omega_trend <- t(chol(omega_trend))
beta_trend <- abs(rcauchy(params$S, 0, scale=0.1))

delta <- marima.sim(kvar=params$S,
                 ar.model=abind(diag(rep(1, params$S)), 
                                diag(beta_trend), along=3),
                 averages=alpha_trend, 
                 resid.cov=crossprod(crossprod(diag(sqrt(theta_trend)), L_omega_trend)),
                 nstart=100, 
                 nsim=periods_to_generate)

delta %>% data.frame() %>%
  set_colnames(paste("series", 1:params$S, sep="_")) %>%
  mutate(period = 1:n()) %>%
  gather(key="series", value="trend", -period) %>%
  ggplot(aes(x = period, y=trend, color=series)) + 
  geom_line(size=0.5, alpha=0.4) +
  geom_vline(data=period_vline, aes(xintercept=xintercept), color="red") +
  theme_minimal() + 
  ggtitle("Generated Trends")
```

## Seasonality

```{r seasonality}
theta_season <- abs(rcauchy(params$S, scale=1)) / params$period_scale
seasoning <- 500
w <- matrix(rnorm(params$S * (periods_to_generate + seasoning), sd=theta_season^2), byrow=T, ncol=params$S)
tau <- c(
  list(
    w[1, ]
  ), 
  list(
    w[2, ] - w[1, ]
  )
)

for (t in 3:(periods_to_generate + seasoning)) {
  if (t < params$s[1]) {
    prior_tau <- do.call(rbind, tau[1:(t-1)])
    tau <- c(tau, list(-colSums(prior_tau) + w[t, ]))
  } else {
    prior_tau <- do.call(rbind, tau[(t - params$s[1] + 1):(t-1)])
    tau <- c(tau, list(-colSums(prior_tau) + w[t, ]))
  }
}

tau <- do.call(rbind, tau)
tau <- tau[(seasoning + 1):nrow(tau), ]

tau %>% data.frame() %>%
  set_colnames(paste("series", 1:params$S, sep="_")) %>%
  mutate(period = 1:n()) %>%
  gather(key="series", value="seasonality", -period) %>%
  ggplot(aes(x = period, y=seasonality, color=series)) + 
  geom_line(size=0.5, alpha=0.4) +
  geom_vline(data=period_vline, aes(xintercept=xintercept), color="red") +
  theme_minimal() + 
  ggtitle("Generated Seasonality")

```

## Cyclicality

```{r cyclicality}
lambda <- 2 * pi / runif(params$S, 250, 1000)
rho <- runif(params$S, 0, 1)
theta_cycle <- abs(rcauchy(params$S, 0, .1)) / params$period_scale
kappa <- matrix(rnorm(params$S * periods_to_generate, sd = theta_cycle), byrow = T, ncol=params$S)
kappa_star <- matrix(rnorm(params$S * periods_to_generate, sd = theta_cycle), byrow = T, ncol=params$S)

omega <- list(
  kappa[1, ]
)
omega_star <- list(
  kappa_star[1, ]
)

for (t in 2:periods_to_generate) { 
  omega <- c(omega, list(
    (rho * cos(lambda) * omega[[t-1]][[1]]) + (rho * sin(lambda) * omega_star[[t-1]][[1]]) + kappa[t, ]
  ))
  omega_star <- c(omega_star, list(
    -(rho * sin(lambda) * omega[[t-1]][[1]]) + (rho * cos(lambda) * omega_star[[t-1]][[1]]) + kappa_star[t, ]
  ))
}

omega <- do.call(rbind, omega)

omega %>% data.frame() %>%
  set_colnames(paste("series", 1:params$S, sep="_")) %>%
  mutate(period = 1:n()) %>%
  gather(key="series", value="cyclicality", -period) %>% 
  ggplot(aes(x = period, y=cyclicality, color=series)) + 
  geom_line(size=0.5, alpha=0.4) +
  geom_vline(data=period_vline, aes(xintercept=xintercept), color="red") +
  theme_minimal() + 
  ggtitle("Generated Cyclicality")
```

## Innovations

```{r innovations}
theta_epsilon <- abs(rcauchy(params$S, 0, .05)) / params$period_scale
omega_epsilon <- rlkjcorr(1, params$S, eta=params$corr_prior)
beta_p <- abs(rnorm(params$S, sd=0.2))
beta_q <- abs(rnorm(params$S, sd=0.2))

garchout <- eccc.sim(periods_to_generate, theta_epsilon, A=diag(beta_p), B=diag(beta_q), 
                     R=omega_epsilon, model="diagonal")

epsilon <- garchout$eps

epsilon %>% data.frame() %>%
  set_colnames(paste("series", 1:params$S, sep="_")) %>%
  mutate(period = 1:n()) %>%
  gather(key="series", value="innovation", -period) %>% 
  ggplot(aes(x = period, y=innovation, color=series)) + 
  geom_vline(data=period_vline, aes(xintercept=xintercept), color="red") +
  geom_line(size=0.5, alpha=0.4) +
  theme_minimal() + 
  ggtitle("Generated Innovations")
```

## View Cumulatives

```{r gen_cumulative}
bind_rows(
  data.frame(delta) %>% mutate(variable="delta"),
  data.frame(tau) %>% mutate(variable="tau"),
  data.frame(omega) %>% mutate(variable="omega"),
  data.frame(epsilon) %>% mutate(variable="epsilon")
) %>%
  set_colnames(c(paste("series", 1:params$S, sep="_"), "variable")) %>%
  group_by(variable) %>% 
  mutate(period = 1:n()) %>%
  gather(key="series", value="trend", -period, -variable) %>%
  group_by(series, variable) %>%
  arrange(period) %>%
  mutate(trend = exp(cumsum(trend))) %>% 
  ggplot(aes(x = period, y=trend, color=series, linetype=variable)) + 
  geom_line(size=0.5, alpha=0.4) +
  geom_vline(data=period_vline, aes(xintercept=xintercept), color="red") +
  theme_minimal() + 
  ggtitle("Generated Trends (Cumulative)")
```

## View Final Price Series

```{r assemble_price_series}
price_fluctuations <- delta + tau + omega + epsilon

starting_prices <- log1p(runif(params$S, 500, 2000))

prices <- rbind(starting_prices, price_fluctuations) %>%
  data.frame() %>%
  set_colnames(paste("series", 1:params$S, sep="_")) %>%
  mutate_all(cumsum) %>%
  dplyr::filter(row_number() != 1) %>% 
  mutate(period = 1:n())

prices %>%
  gather(key="series", value="log_price", -period) %>%
  mutate(price = expm1(log_price)) %>%
  ggplot(aes(x = period, y = price, color=series)) + 
  geom_line(size=0.5, alpha=0.7) + 
  geom_vline(data=period_vline, aes(xintercept=xintercept), color="red") +
  theme_minimal() + 
  ggtitle("Generated Price Series")
```


# Make Model And Simulate

```{r assemble_stan_data_function}
assemble_stan_data <- function(prices, predictors) compose_data(
                          # The data
                          prices, 
                          x=predictors, 
                          N_periods = max(prices$period),
                          N_features = ncol(predictors),
                          # Parameters controlling the model 
                          periods_to_predict = params$periods_to_predict, 
                          p = params$p, 
                          q = params$q, 
                          ar = params$ar, 
                          period_scale = params$period_scale, 
                          s = array(params$s, dim = length(params$s)),
                          N_seasonality = length(params$s),
                          cyclicality_prior = 52 * 5, 
                          corr_prior = params$corr_prior,
                          # Parameters controlling hierarchical shrinkage
                          nu = params$nu, 
                          m0 = params$m0,
                          .n_name = n_prefix("N"))
```

```{r make_stan_data}
prices_gathered <- prices %>%
  gather(key="series", value="log_price", -period) %>%
  mutate(price = exp(log_price) - 1) %>%
  dplyr::select(-log_price) %>%
  arrange(period, series) %>%
  mutate(id = 1:n())

sampled_prices_for_stan <- data.frame(
  period = sample(params$P, params$N, replace=T), 
  series = sample(unique(prices_gathered$series), params$N, replace=T)
) %>% 
  inner_join(prices_gathered) %>%
  mutate(
    y = rnorm(n(), mean=price, sd=50)
  ) %>%
  dplyr::filter(y > 0) %>%
  mutate(
    weight = 1,
    id = 1:n()
  )

predictors <- matrix(rnorm(params$P * params$F), ncol=params$F)

sampled_stan_data <- assemble_stan_data(dplyr::select(sampled_prices_for_stan, -id), predictors)

exact_stan_data <- prices_gathered %>% 
  dplyr::filter(period <= params$P) %>% 
  rename(y=price) %>% 
  mutate(weight = 1) %>% 
  assemble_stan_data(predictors)

prices_gathered %<>% rename(y=price)
```

```{r select_price_set}
exact <- ! as.logical(params$use_sampled)

if (exact) {
  stan_data <<- exact_stan_data
  price_set <<- prices_gathered
} else {
  stan_data <<- sampled_stan_data
  price_set <<- sampled_prices_for_stan 
}
```


```{r make_model}
model <- stan_model("./mbsts.stan")
```

```{r sample}
init_f <- function(chain_id=1) {
  list(alpha_ar="0")
}
samples <- sampling(model, data=stan_data, 
                    chains=3, 
                    #warmup=1000,
                    #iter=3000, thin=3, init=init_f, 
                    cores=5, save_warmup=FALSE,
                    refresh=200,
                    ) %>% recover_types(price_set)
```

# Inspect Model

## Diagnostics

```{r check_hmc}
check_hmc_diagnostics(samples)
```

```{r energy}
bayesplot::mcmc_nuts_energy(nuts_params(samples))
```

## Pairs

```{r pairs}
parameters <- samples %>% gather_draws(
  alpha_ar[series], 
  beta_ar[feature, series], 
  delta_t0[feature, series], 
  theta_ar[series], 
  theta_season[series, seasonality], 
  lambda[series], 
  rho[series], 
  theta_cycle[series], 
  beta_xi[feature, series], 
  omega_garch[series], 
  beta_p[feature, series], 
  beta_q[feature, series], 
  starting_prices[series]
)
```

```{r pairs,eval=F}
series_to_show <- "series_1"

parameters %>%
  dplyr::filter(series == "series_1") %>% 
  unite("variable", .variable, feature, seasonality) %>%
  mutate(
    variable = str_replace_all(variable, "_NA", "")
  ) %>%
  spread(key="variable", value=".value") %>% 
  ungroup() %>% 
  dplyr::select(-series, -.chain, -.iteration, -.draw) %>%
  ggpairs(
    upper=list(continuous="density"), 
    lower=list(continuous="smooth_loess"), 
    diag=list(continuous="densityDiag"), 
    progress=FALSE
  )
```

## Prediction Accuracy

### Residuals 

```{r recover_y_hat}
y_hat <- samples %>% gather_draws(log_y_hat[id]) %>%
  mutate(.value = expm1(.value))
```

```{r residuals}
y_hat %>% 
  dplyr::select(id, .value) %>%
  inner_join(price_set) %>% 
  mutate(
    error = .value - y
  ) %>%
  group_by(id, series) %>%
  summarize(
    mean_error = mean(error),
    rmse = sqrt(mean(error^2)), 
    sd_y_hat = sd(.value)
  ) %>% 
  dplyr::select(series, mean_error, rmse, sd_y_hat) %>%
  gather(key="metric", value="value", -series, -id) %>% 
  ggplot(aes(x = value, fill=series)) +
  geom_histogram(bins=40, position="identity", alpha = 1/params$S) + 
  facet_grid(series ~ metric, scales="free") +
  theme_minimal() +
  ggtitle("Residuals")
```

```{r visualize_predictions}
price_set %>% dplyr::select(id, y, series) %>%
  inner_join(y_hat %>% 
               group_by(id) %>%
               summarize(y_hat = mean(.value))) %>%
  ggplot(aes(x = y, y=y_hat)) +
  geom_point(size=0.2, alpha=0.2) +
  facet_wrap(~ series) +
  scale_x_continuous(limits=c(0, NA)) + 
  scale_y_continuous(limits=c(0, NA)) + 
  annotate("segment", x = 0, y = 0, xend=max(price_set$y), yend=max(price_set$y), color="red", size=0.5, alpha=0.5) +
  theme_minimal() + 
  ggtitle("Prices vs. Mean Predictions")
```

### Recovered Time Series

```{r recover_time_series}
periodic_price_hat <- samples %>% gather_draws(log_prices[period, series])
```

```{r time_series_by_chain}
mean_periodic_price_hat <- periodic_price_hat %>%
  dplyr::filter(is.finite(expm1(.value))) %>% 
  group_by(period, series) %>%
  summarize(
    y_hat = expm1(mean(.value))
  ) %>%
  mutate(.chain = as.character("mean_hat"))

periodic_price_hat %>%
  dplyr::filter(is.finite(.value)) %>% 
  group_by(period, series, .chain) %>%
  summarize(
    y_hat = expm1(mean(.value))
  ) %>%
  ungroup() %>%
  mutate(.chain = as.character(.chain)) %>%
  bind_rows(prices_gathered %>% mutate(.chain="gen") %>% rename(y_hat=y), mean_periodic_price_hat) %>% 
  mutate(
#    alpha = case_when(
#      .chain %in% c("gen", "mean_hat") ~ .8, 
#      TRUE ~ .2
#    ), 
    .chain = fct_relevel(factor(.chain), c("gen", "mean_hat"))
  ) %>% 
  ggplot(aes(x = period, y = y_hat, color=series, linetype=.chain)) +
  geom_line(size=0.3, alpha=0.6) +
  theme_minimal() +
  ggtitle("Reconstruction of Time Series By Chain")
```

### Forecast

```{r gather_price_forecast}
price_forecast <- samples %>% 
  gather_draws(log_prices_hat[period, series]) 
```

```{r visualize_price_forceast}
true_future <- prices_gathered %>%
  dplyr::filter(period > params$P) %>%
  mutate(period = period - params$P)
  
price_forecast %>%
  dplyr::filter(is.finite(expm1(.value))) %>%
  ungroup() %>%
  group_by(series, period, .chain) %>%
  summarize(y = expm1(mean(.value))) %>% 
  #mutate(chain = "mean_hat") %>%
  mutate(.chain = as.character(.chain)) %>%
  bind_rows(true_future %>% dplyr::select(-id) %>% mutate(.chain="gen")) %>%
  ggplot(aes(x=period, y=y, color=series, linetype=.chain)) +
  geom_line() + 
  geom_lineribbon(alpha=0.3, size=0.5, data=price_forecast %>% 
                    dplyr::filter(is.finite(.value)) %>%
                    mutate(y=expm1(.value), .chain="median") %>%
                    ungroup() %>%
                    dplyr::select(series, period, y, .chain) %>% 
                    group_by(period, series, .chain) %>% 
                    median_qi(y, .width=c(0.05, 0.1, 0.2))
                    ) +
  scale_fill_brewer("Confidence", palette="Blues") +
  theme_minimal() +
  scale_y_log10() +
  facet_grid(series ~ ., scales="free") +
  ggtitle("Forecasted Prices")
```

### Combined

```{r combined_plot}
bind_rows(prices_gathered %>% 
            mutate(.chain="Actual") %>% 
            rename(y_hat=y), 
          mean_periodic_price_hat) %>% 
  #dplyr::filter(.chain %in% c("forecast")) %>% 
  ggplot(aes(x = period, y = y_hat, color=series, linetype=.chain)) +
  geom_lineribbon(alpha=0.2, size=0.5, data=price_forecast %>% 
                  dplyr::filter(is.finite(.value)) %>%
                  ungroup() %>%
                  mutate(
                    y_hat=expm1(.value), 
                    .chain="median",
                    period = period + max(mean_periodic_price_hat$period)
                  ) %>%
                  dplyr::select(series, period, y_hat, .chain) %>% 
                  group_by(period, series, .chain) %>% 
                  median_qi(y_hat, .width=c(0.05, 0.1, 0.2))
                  ) +
  geom_line(size=0.5, alpha=0.8) +
  scale_fill_brewer("Confidence", palette="Blues") +
  theme_minimal() +
  ggtitle("Reconstructed Time Series and Forecast vs. Actual", 
          subtitle=glue("AR=", stan_data$ar, ", P=", stan_data$p, ", Q=", stan_data$q))
```

## Components 

### Visualize All Components

```{r component_view_functions}
components_by_mean <- function(components) {
  components %>%
  ungroup() %>%
  group_by(period, series, component) %>%
  summarize(.value = mean(.value)) %>% 
  mutate(.chain = "mean_hat") %>%
  ggplot(aes(x = period, y = .value, color=component)) +
  facet_grid(series ~ ., scales="free_y") +
  geom_line(size=0.2) +
  theme_minimal() +
  ggtitle("Components of Price Movements (Means)")
}

components_by_chain <- function(components) {
  components %>%
    ungroup() %>%
    group_by(period, series, component, .chain) %>%
    summarize(.value = mean(.value)) %>%
    mutate(.chain = as.character(.chain)) %>%
    bind_rows(
      components %>%
        ungroup() %>%
        group_by(period, series, component) %>%
        summarize(.value = mean(.value)) %>% 
        mutate(.chain = "mean_hat")
    ) %>%
    ggplot(aes(x = period, y = .value, color=component, linetype=factor(.chain))) +
    facet_grid(series ~ ., scales="free_y") +
    geom_line(size=0.2) +
    theme_minimal() +
    ggtitle("Components of Price Movements (By Chain)")
}

cumulative_components <- function(components) {
  components %>%
    ungroup() %>%
    group_by(period, series, component, .chain) %>%
    summarize(.value = mean(.value)) %>%
    mutate(.chain = as.character(.chain)) %>%
    bind_rows(
      components %>%
        ungroup() %>%
        group_by(period, series, component) %>%
        summarize(.value = mean(.value)) %>% 
        mutate(.chain = "mean_hat")
    ) %>%
    group_by(component, .chain, series) %>%
    arrange(period) %>%
    mutate(.value = cumsum(.value)) %>% 
    ggplot(aes(x = period, y = .value, color=component, linetype=factor(.chain))) +
    facet_grid(series ~ ., scales="free_y") +
    geom_line(size=0.2) +
    #scale_x_continuous(limits=c(0, 80)) +
    theme_minimal() +
    ggtitle("Cumulative Components of Price Movements (By Chain)")
}

```

```{r recover_components}
combined_components <- samples %>%  gather_draws(delta[period, series], tau_s[seasonality][period, series], omega[period, series], epsilon[period, series], xi[period, series])  %>% 
  ungroup() %>% 
  mutate(
    component = case_when(
      .variable == "delta" ~ "trends", 
      .variable == "epsilon" ~ "innovations", 
      .variable == "tau_s" ~ paste("seasonality", seasonality, sep="_"), 
      .variable == "omega" ~ "cyclicality", 
      .variable == "xi" ~ "predictors"
    )
  )

components_by_mean(combined_components)
```

```{r view_components_by_chain}
components_by_chain(combined_components)
```

```{r components_cumulative}
cumulative_components(combined_components)
```

## Forecast Components

```{r recover_fcomponents}
forecast_components <- samples %>%  gather_draws(delta_hat[period, series], tau_hat_all[period, series], omega_hat[period, series], epsilon_hat[period, series])  %>% 
  ungroup() %>% 
  mutate(
    component = case_when(
      .variable == "delta_hat" ~ "trends", 
      .variable == "epsilon_hat" ~ "innovations", 
      .variable == "tau_hat_all" ~ "seasonality",
      .variable == "omega_hat" ~ "cyclicality", 
      .variable == "xi_hat" ~ "predictors"
    )
  )

components_by_mean(forecast_components)
```


```{r view_fcomponents_by_chain}
components_by_chain(forecast_components)
```


```{r fcomponents_cumulative}
cumulative_components(forecast_components)
```

## Components in Detail

### Cyclicality

```{r cyclicality}
cylicality_hat <- samples %>% gather_draws(lambda[series], rho[series])
```

```{r view_cyclicality}
true_cyclicality <- data.frame(
  lambda = lambda, 
  rho = rho, 
  series = sort(unique(cylicality_hat$series))
) %>%
  gather(key=".variable", value=".value", -series)

cylicality_hat %>% 
  ggplot(aes(x = .value)) + 
  geom_histogram(bins=60) +
  geom_vline(data=true_cyclicality, aes(xintercept=.value), color="blue", size=0.5) +
  facet_grid(series ~ .variable, scales="free") +
  scale_x_log10() +
  theme_minimal() + 
  ggtitle("Estimates of Cycle Frequency and Damping")
```

## Correlations 

### Trend Correlations

```{r recover_return_correlations}
trend_corr_hat <- samples %>% gather_draws(trend_corr[series1, series2])

true_trends <- omega_trend %>%
  data.frame() %>%
  set_colnames(1:ncol(omega_trend)) %>%
  mutate(series1 = 1:n()) %>%
  gather(key="series2", value=".value", -series1)
```

```{r show_trend_corr}
trend_corr_hat %>%
  dplyr::filter(series1 > series2) %>%
  ggplot(aes(x = .value)) +
  geom_histogram(bins=100) +
  geom_vline(data=data.frame(xintercept=0), aes(xintercept=xintercept), color="red", size=0.2) +
  geom_vline(data=true_trends, aes(xintercept=.value), color="blue", size=0.2) +
  facet_grid(series1 ~ series2, scales="free_y") +
  theme_minimal() + 
  scale_x_continuous(limits=c(-1, 1)) +
  ggtitle("Trend Correlations")
```

### Innovation Correlations

```{r recover_innovation_correlations}
innovation_corr_hat <- samples %>% gather_draws(innovation_corr[series1, series2])

true_innovations <- omega_epsilon %>%
  data.frame() %>%
  set_colnames(1:ncol(omega_trend)) %>%
  mutate(series1 = 1:n()) %>%
  gather(key="series2", value=".value", -series1)
```

```{r show_innovation_corr}
innovation_corr_hat %>%
  dplyr::filter(series1 > series2) %>%
  ggplot(aes(x = .value)) +
  geom_histogram(bins=100) +
  geom_vline(data=data.frame(xintercept=0), aes(xintercept=xintercept), color="red", size=0.2) +
  geom_vline(data=true_innovations, aes(xintercept=.value), color="blue", size=0.2) +
  facet_grid(series1 ~ series2) +
  theme_minimal() + 
  ggtitle("Innovation Correlations")
```


## Coefficients and Shrinkage Effects

### Visualize Coefficients

```{r visualize_coefficients}
coefficients_hat <- samples %>% gather_draws(beta_ar_c[feature, series], beta_xi[feature, series], 
                                             beta_p_c[feature, series], beta_q_c[feature, series])

coefficients_hat %>%
  ggplot(aes(x = .value)) +
  geom_histogram(bins=60) +
  facet_grid(series ~ .variable + feature, scales="free_x") +
  theme_minimal() +
  ggtitle("Estimated Coefficients")
```


## Volatilities and Intercepts 

### Visualize Volatilities and Intercepts

```{r visualize_volatilities}
volatilities_hat <- samples %>% gather_draws(theta_ar[series], alpha_ar[series], theta_cycle[series], theta_season[seasonality][series], omega_garch[series]) %>%
  ungroup() %>% 
  mutate(
    component = case_when(
      .variable == "theta_season" ~ paste("theta_season", seasonality, sep="_"),
      TRUE ~ .variable
    )
  )

volatilities_hat %>% 
  ggplot(aes(x = .value)) +
  geom_histogram(bins=60) +
  facet_grid(component ~ series, scales="free_x") +
  theme_minimal() +
  ggtitle("Estimated Volatilities and Baselines")
```


## Further Inspection

```{r shinystan}
library(shinystan)
ss_samples <- samples %>% as.shinystan() %>%
  drop_parameters(pars=c("log_prices_hat", "delta_hat", "tau_hat_all", 
                  "omega_hat", "omega_star_hat", "theta_hat", 
                  "epsilon_hat", "nu_ar_hat", "kappa_hat", 
                  "kappa_star_hat", "w_t_hat", 
                  "log_prices", "delta", "tau_s", "tau", "omega", "omega_star", 
                  "theta", "xi", 
                  "nu_trend", "w_t", "kappa", "kappa_star", "epsilon"))
```


# DONE