time series

timeseries.R

# @ericpgreen

# CausalImpact
  library(CausalImpact)
  set.seed(1)
  x1 <- 24 + arima.sim(model = list(ar = 0.999), n = 24)
  y <- 1.2 * x1 + rnorm(24)
  y[13:24] <- y[13:24] + 10
  data <- cbind(y, x1)
  
  matplot(data, type = "l")
  
  pre.period <- c(1, 12)
  post.period <- c(13, 24)
  
  impact <- CausalImpact(data, pre.period, post.period)
  plot(impact)
  
  time.points <- seq.Date(as.Date("2018-01-01"), by = "month", length.out = 24)
  data <- zoo(cbind(y, x1), time.points)
  head(data)

  summary(impact)
  
  # Posterior inference {CausalImpact}
  # 
  # Average      Cumulative  
  # Actual                   22           268         
  # Prediction (s.d.)        12 (0.35)    143 (4.17)  
  # 95% CI                   [11, 13]     [134, 151]  
  # 
  # Absolute effect (s.d.)   10 (0.35)    126 (4.17)  
  # 95% CI                   [9.8, 11]    [117.4, 134]
  # 
  # Relative effect (s.d.)   88% (2.9%)   88% (2.9%)  
  # 95% CI                   [82%, 94%]   [82%, 94%]  
  # 
  # Posterior tail-area probability p:   0.001
  # Posterior prob. of a causal effect:  99.8999%
  # 
  # For more details, type: summary(impact, "report")
  
  impact$summary
  
  #               Actual      Pred  Pred.lower  Pred.upper    Pred.sd AbsEffect
  # Average     22.34431  11.87668    11.17645    12.56069  0.3474998  10.46763
  # Cumulative 268.13172 142.52016  134.11738  150.72834 4.1699972 125.61156
  

# its.analysis
  library(its.analysis)
  library(tidyverse)
  data2 <- as.data.frame(data)
  data2 <- tibble::rownames_to_column(data2, "time")
  data2$int <- c(rep(0, 12), rep(1, 12))
  
  itsa.model(data=data2, time="time", depvar="y", interrupt_var = "int", 
             covariates = "x1", alpha=0.05, bootstrap=TRUE, Reps = 250)


  # $group.means
  # interrupt_var count     mean     s.d.
  # 1             0    12 11.68917 1.848760
  # 2             1    12 22.34431 2.666517
  
  # is it true that we compare results as:
  
  # CausalImpact / its.analysis
  # 0 11.88 / 11.69
  # 1 22.34 / 22.34