library("CausalImpact")
library("tidyverse")
Over the next few weeks, I want to learn more about Bayesian Neural Networks. Here in this blog post, I will gather resources and examples.
Causal Impact
I have seen several people mention (on Twitter or Slack) the CausalImpact
package. Here, I roughly follow the vignette.
Toy Data
set.seed(320)
# create time series and apply vertical shift
<- 100 + arima.sim(model = list(ar = 0.999), n = 100)
ts_raw
# add "noise"
<- 1.2*ts_raw + rnorm(100)
y
# intentionlly "break" time series and lift a part
71:100] <- y[71:100] + 25
y[
<- 1:100
x <- data.frame(x,y) df
|>
df ggplot(aes(x,y)) +
geom_vline(xintercept = 70, color = "gray50",
linetype = 2, linewidth = 3) +
geom_line(color = "red", linewidth = 2) +
labs(title = "Time Series",
subtitle = "with an induced discontinuity") +
theme_minimal()
Bayesian Analysis
training over pre-intervention period
testing over post-intervention period
the
CausalImpact
function- assembles structural time-series model
- performs posterior inference
- computes estimates for casual effect
<- c(1,70)
pre_intervention <- c(71, 100)
post_intervention <- CausalImpact::CausalImpact(df,
ts_impact
pre_intervention, post_intervention)
Model Statistics
summary(ts_impact)
Posterior inference {CausalImpact}
Average Cumulative
Actual 86 2565
Prediction (s.d.) 103 (3.8) 3100 (115.3)
95% CI [96, 111] [2882, 3326]
Absolute effect (s.d.) -18 (3.8) -535 (115.3)
95% CI [-25, -11] [-761, -317]
Relative effect (s.d.) -17% (3.1%) -17% (3.1%)
95% CI [-23%, -11%] [-23%, -11%]
Posterior tail-area probability p: 0.00106
Posterior prob. of a causal effect: 99.89418%
For more details, type: summary(impact, "report")
Visualization
The plot
of an “impact” object returns a ggplot
object!
plot(ts_impact)
bsts
Bayesian structural time series
- looking at the first example on this blog post
library("bsts")
data(iclaims) #unemployment data
Trend and Seasonal Components
<-
state_specs ::AddSemilocalLinearTrend(list(),
bsts$iclaimsNSA)
initial.claims<- bsts::AddSeasonal(state_specs,
state_specs $iclaimsNSA,
initial.claimsnseasons = 52)
Model
<- bsts::bsts(initial.claims$iclaimsNSA,
model_1 state.specification = state_specs,
niter = 100)
=-=-=-=-= Iteration 0 Tue Mar 19 23:21:31 2024
=-=-=-=-=
=-=-=-=-= Iteration 10 Tue Mar 19 23:21:31 2024
=-=-=-=-=
=-=-=-=-= Iteration 20 Tue Mar 19 23:21:31 2024
=-=-=-=-=
=-=-=-=-= Iteration 30 Tue Mar 19 23:21:32 2024
=-=-=-=-=
=-=-=-=-= Iteration 40 Tue Mar 19 23:21:32 2024
=-=-=-=-=
=-=-=-=-= Iteration 50 Tue Mar 19 23:21:32 2024
=-=-=-=-=
=-=-=-=-= Iteration 60 Tue Mar 19 23:21:32 2024
=-=-=-=-=
=-=-=-=-= Iteration 70 Tue Mar 19 23:21:33 2024
=-=-=-=-=
=-=-=-=-= Iteration 80 Tue Mar 19 23:21:33 2024
=-=-=-=-=
=-=-=-=-= Iteration 90 Tue Mar 19 23:21:33 2024
=-=-=-=-=
Visualization
plot(model_1)
Prediction
- predict next 12 time points
- along last 156 time points (3 years)
<- predict(model_1, horizon = 12)
pred_3y plot(pred_3y, plot.original = 156)
Online Videos
- Bayesian Neural Network by TwinEd Productions (7 minute video)