dataTS$Y <- round(dataTS$Y,1)

d.temp <- rbind( head(dataTS), c("...","...","...","..."),
                 dataTS[ 198:200, ], 
                 c("**START**","**TREATMENT**","---","---"), 
                 dataTS[ 201:204, ],  
                 c("...","...","...","..."),
                 tail(dataTS) )

row.names( d.temp ) <- NULL  
pander( d.temp )
Y T D P
38.5 1 0 0
48.6 2 0 0
102.4 3 0 0
58.1 4 0 0
60 5 0 0
107.8 6 0 0
… … … …
60 198 0 0
79.4 199 0 0
62.4 200 0 0
START TREATMENT — —
173.1 201 1 1
147.3 202 1 2
100.7 203 1 3
125.7 204 1 4
… … … …
300 360 1 160
217.5 361 1 161
243.9 362 1 162
233.3 363 1 163
256.6 364 1 164
251.1 365 1 165 
plot( dataTS$T, dataTS$Y,
      bty="n", pch=19, col="gray",
      ylim = c(0, 300), xlim=c(0,400),
      xlab = "Time (days)", 
      ylab = "Wellbeing Index" )

# Line marking the interruption
abline( v=200, col="firebrick", lty=2 )
text( 200, 300, "Start of Wellbeing Classes", col="firebrick", cex=1.3, pos=4 )

# Add the regression line
ts <- lm( Y ~ T + D + P, data=dataTS )
lines( dataTS$T, ts$fitted.values, col="steelblue", lwd=2 )

regTS = lm ( Y ~ T + D + P, data=dataTS)  # Our time series model
stargazer( regTS, 
           type = "html", 
           dep.var.labels = ("Wellbeing"),
           column.labels = ("Model results"),
           covariate.labels = c("Time", "Treatment", 
                                "Time Since Treatment"),
           omit.stat = "all", 
           digits = 2 )
Dependent variable:
Wellbeing
Model results
Time 0.19***
(0.04)
Treatment 13.10**
(6.12)
Time Since Treatment 0.54***
(0.06)
Constant 57.15***
(4.13)
Note: p<0.1; p<0.05; p<0.01 
# We create a small dataset with the new values
data1 <- as.data.frame( cbind( T = 201, D = 1, P = 1 )) 

# We use the function predict to (1) take the 
#  coefficients estimated in regTS and 
#  (2) calculate the outcome Y based on the 
#  values we set in the new datset
y1 <- predict( regTS, data1 ) 

# We plot our initial observations, the column Y in our dataset
plot( dataTS$Y,
      bty="n",
      col = gray(0.5,0.5), pch=19,
      xlim = c(1, 365), 
      ylim = c(0, 300),
      xlab = "Time (days)", 
      ylab = "Wellbeing Index")

# We add a point showing the level of wellbeing at time = 201)
points( 201, y1, col = "dodgerblue4", 
        pch = 19, bg = "dodgerblue4", cex = 2 )
text( 201, y1, labels = "t = 201", pos = 4, cex = 1 )

# Line marking the interruption
abline( v=200, col="red", lty=2 )

# Figure 2.2: Wellbeing level at t=230
# Add the new point at t=230 to the graph
data2 <- as.data.frame( cbind( T = 230, D = 1, P = 30 )) # New data

y2 <- predict( regTS, data2 ) # We predict the new outcome

# We plot our initial observations, the column Y in our dataset
plot( dataTS$Y,
      bty="n",
      col = gray(0.5,0.5), pch=19,
      xlim = c(1, 365), 
      ylim = c(0, 300),
      xlab = "Time (days)", 
      ylab = "Wellbeing Index")

# We add a point showing the level of wellbeing at time = 201
points(201, y1, col = "dodgerblue4", pch = 19, bg = "red", cex = 2)

# We add a new point showing the level of wellbeing at time = 230)
points(230, y2, col = "dodgerblue4", pch = 19, bg = "red", cex = 2)

# Label for our predicted outcome
text(201, y1, labels = "t = 201", pos = 4, cex = 1)

#Label for the counterfactual 
text(230, y2, labels = "t = 230", pos = 4, cex = 1)

# Line marking the interruption
abline( v=200, col="red", lty=2 )

# The Counterfactual
# NOte: Counterfactual varies for each point in time. Plot a new point at time = 320 and its counterfactual
data3 <- as.data.frame(cbind( T= 230, D = 0, P = 0)) # Data if the intervention does not occur

y3 <- predict(regTS, data3) #Counterfactual

# We plot our initial observations, the column Y in our dataset
plot( dataTS$Y,
      bty="n",
      col = gray(0.5,0.5), pch=19,
      xlim = c(1, 365), 
      ylim = c(0, 300),
      xlab = "Time (days)", 
      ylab = "Wellbeing Index")

# We add a  point showing the level of wellbeing at time = 230
points(230, y2, col = "dodgerblue4", pch = 19, bg = "red", cex = 2)

# We add a point indicating the counterfactual
points(230, y3, col = "darkorange2", pch = 19, cex = 2)

# Label for our predicted outcome
text(230, y2, labels = "Y at t = 230", pos = 4, cex = 1)

#Label for the counterfactual 
text(230, y3, labels = "C at t = 230", pos = 4, cex = 1)

# Line marking the interruption
abline( v=200, col="red", lty=2 )

# Counterfactual varies for each point in time. Plot a new point at time = 320 and its counterfactual
data4 <- as.data.frame(cbind( T = 320, D = 1, P = 80)) 
data5 <- as.data.frame(cbind( T = 320, D = 0, P = 0))

y4 <- predict(regTS, data4)
y5 <- predict(regTS, data5)

# We plot our initial observations, the column Y in our dataset
plot( dataTS$Y,
      bty="n",
      col = gray(0.5,0.5), pch=19,
      xlim = c(1, 365), 
      ylim = c(0, 300),
      xlab = "Time (days)", 
      ylab = "Wellbeing index")

# Wellbeing at time = 230
points(230, y2, col = "dodgerblue4", pch = 19, bg = "red", cex = 2)

# Counterfactual at time = 230
points(230, y3, col = "darkorange2", pch = 19, cex = 2)

# Wellbeing at time = 320
points(320, y4, col = "dodgerblue4", pch = 19, cex = 2)

# Counterfactual at time = 320
points(320, y5, col = "darkorange2", pch = 19, cex = 2)

# Adding labels
text(320, y4, labels = "Y at t = 320", pos = 4, cex = 1)
text(320, y5, labels = "C at t = 320", pos = 4, cex = 1)
text(230, y2, labels = "Y at at = 230", pos = 4, cex = 1)
text(230, y3, labels = "C at t = 230", pos = 4, cex = 1)

# Line marking the interruption
abline( v=200, col="red", lty=2 )

# Finally, we can look at the results by plotting all predicted outcomes and their counterfactuals
pred1 <- predict(regTS, dataTS) 
# To estimate all predicted values of Y, we just use our dataset

datanew <- as.data.frame(cbind(T = rep(1 : 365), D = rep(0), P = rep(0))) 
# Create a new dataset where Treatment and Time Since Treatment are equal to 0 as the intervention did not occur.

pred2 <- predict(regTS, datanew) 
# Predict the counterfactuals

plot( dataTS$Y,
      bty="n",
      col = gray(0.5,0.5), pch=19,
      xlim = c(1, 365), 
      ylim = c(0, 400),
      xlab = "Time (days)", 
      ylab = "Wellbeing index")

lines( rep(1:199), pred1[1:199], col="dodgerblue4", lwd = 3 )
lines( rep(201:365), pred1[201:365], col="dodgerblue4", lwd = 3 )
lines( rep(200:365), pred2[200:365], col="darkorange2", lwd = 3, lty = 5 ) 

text(0, 45, labels = "Predicted values", pos = 4, cex = 1, col = "dodgerblue3")
text(300, 95, labels = "Counterfactual", pos = 4, cex = 1, col = "darkorange2")

# Line marking the interruption
abline( v=200, col="darkorange2", lty=2 )

# Delayed effects in time series
# Note: The model coefficients above does not tell us if the difference between predicted outcomes and its counterfactuals is statisticaaly significant
# It only tells us:
# 1. There is an immediate change after the intervention and
# 2. The slope has changed after the intervention
reg2 = lm(Y ~ T + D + P, data = dataTS)
stargazer( reg2, 
           type = "html", 
           dep.var.labels = ("Wellbeing"),
           column.labels = ("Model results"),
           covariate.labels = c("Time", "Treatment", "Time Since Treatment"),
           omit.stat = "all", 
           digits = 2 )
Dependent variable:
Wellbeing
Model results
Time 0.19***
(0.04)
Treatment 13.10**
(6.12)
Time Since Treatment 0.54***
(0.06)
Constant 57.15***
(4.13)
Note: p<0.1; p<0.05; p<0.01 
# Delayed time effect presented in a graph
# Thinking about when the effect is likely to occur is important in time series
# Figure: The effect of the policy changes over time
pred1 <- predict(reg2, dataTS) 

datanew <- as.data.frame(cbind(T = rep(1 : 365), D = rep(0), P = rep(0))) 

pred2 <- predict(reg2, datanew) 

plot(dataTS$Y,
    col = "gray",
    xlim = c(1, 365), 
    ylim = c(0, 400),
    xlab = "Time (days)", 
    ylab = "Wellbeing index")

lines( rep(1:199), pred1[1:199], col="dodgerblue4", lwd = 3 )
lines( rep(201:365), pred1[201:365], col="dodgerblue4", lwd = 3 )
lines( rep(200:365), pred2[200:365], col="darkorange2", lwd = 3, lty = 5 ) 

text(0, 45, labels = "Predicted values", pos = 4, cex = 1, col = "dodgerblue3")
text(300, 105, labels = "Counterfactual", pos = 4, cex = 1, col = "darkorange2")

# Line marking the interruption
abline( v=200, col="red", lty=2 )

# Line at t = 300
abline( v=300, col="forestgreen", lty=2 )