rm(list=ls()) #clear the workspace
set.seed(2020) #so that everyone gets the same results

#Set the working directory to the folder containing W3_BLIK_SCM.csv, e.g.:
#In RStudio: Session -> Set Working Directory -> To Source File Location

#note: the script uses the pipe operator |> (requires R 4.1 or newer)
#note: require() alone does not load the package after installing it - hence library() in the braces
if(!require("rddtools")) {install.packages("rddtools"); library("rddtools")}
if(!require("tidysynth")) {install.packages("tidysynth"); library("tidysynth")}


#####################################
#1. Difference-in-differences - illustration
#we draw the plot for the control group
plot(c(0, 1), c(5, 7), 
     type = "p",
     ylim = c(4, 12),
     xlim = c(-0.2, 1.2),
     main = "DID estimator",
     xlab = "Time",
     ylab = "y",
     col = "grey",
     pch = 19,  
     cex = 1.2,
     xaxt = "n",
     yaxt = "n")


axis(1, at = c(0, 1), labels = c("before", "after"))


#the treated group
points(c(0, 1), c(7, 11), 
       col = "blue",
       pch = 19,
       cex = 1.2)

#the treated group had it followed the same trend as the control group
points(1,9, 
       col = "lightblue",
       pch = 19,
       cex = 1.2)

#connect the points
lines(c(0, 1), c(7, 11), col = "blue")
lines(c(0, 1), c(5, 7), col = "grey")
lines(c(0, 1), c(7, 9), col = "lightblue", lty = 2)
lines(c(1, 1), c(9, 11), col = "black", lty = 2, lwd = 2)

#add the annotation
text(1, 9.7, expression(hat(beta)[1]^{DID}), cex = 0.8, pos = 4)
text(1, 10.3, "treatment effect", cex = 0.8, pos = 4)

#now an example on generated data
N<-250 #set the sample size

#define the treatment effect
TreatEff<-5

#create a variable that will split the N observations into groups
S<- c(rep(0, N/2), rep(1, N/2))

#create the before- and after-treatment data
#note the parentheses: 1:(N/2) is observations 1-125, while 1:N/2 is (1:N)/2 - a completely different thing!
y_przed <- 7 + rnorm(N)
y_przed[1:(N/2)] <- y_przed[1:(N/2)] - 1
y_po <- 7 + 2 + TreatEff * S + rnorm(N)
y_po[1:(N/2)] <- y_po[1:(N/2)] - 1 

     
#let's compute the effect
mean(y_po[S == 1]) - mean(y_przed[S == 1]) - (mean(y_po[S == 0]) - mean(y_przed[S == 0]))
#why didn't we get exactly 5?

#let's compute it using OLS
summary(lm(y_po - y_przed ~ S))


#####################################
#2. DiD with staggered treatment timing (staggered adoption)
#We show why a plain TWFE regression can fail when the intervention takes
#effect at different times and the effects are heterogeneous.

N_j <- 300; T_max <- 20
panel <- expand.grid(id = 1:N_j, t = 1:T_max)

#three cohorts are treated at different moments
kohorta <- sample(c(5, 10, 15), N_j, replace = TRUE)
panel$g <- kohorta[panel$id]

#the TRUE effect: grows over time since treatment and differs across cohorts
#(the earlier the cohort, the stronger the effect - this alone is enough to distort TWFE)
panel$staz <- pmax(0, panel$t - panel$g)
sila <- c("5" = 3.0, "10" = 1.5, "15" = 0.5)
panel$efekt <- sila[as.character(panel$g)] * panel$staz
panel$D <- as.numeric(panel$t >= panel$g)

panel$y <- panel$id*0.01 + panel$t*0.2 + panel$efekt + rnorm(nrow(panel))

#true average effect among the treated observations
prawdziwy <- mean(panel$efekt[panel$D == 1])

#the TWFE estimator (unit and time fixed effects)
twfe <- lm(y ~ D + factor(id) + factor(t), data = panel)
oszacowany <- coef(twfe)["D"]

c(true = prawdziwy, TWFE = oszacowany)
#TWFE does not recover the true effect - some of the comparisons use units that are
#ALREADY treated as the control group (Goodman-Bacon 2021)

#Solutions: the did package (Callaway, Sant'Anna), fixest::sunab (Sun, Abraham),
#didimputation. See also the event study plot: fixest::feols + iplot


#####################################
#3. Regression discontinuity design
#example: university admission decided by a test score (70% threshold)
x <- runif(1000, 0, 1)
y <- as.numeric(x >= 0.7) #equivalent to a loop: 1 if the threshold is crossed, 0 otherwise

plot(x,y,
     col = "blue",
     cex = 0.35,
     xlab = "Test score", 
     ylab = "University admission")

#Technical note: the rdd package (the RDestimate function) was removed from CRAN on 10 July 2025.
#We use the rddtools package, which returned to CRAN on 29 October 2025 (version 2.0.2)
#and offers BOTH non-parametric AND parametric estimation.
#An alternative used in the literature today: the rdrobust package (Calonico, Cattaneo, Titiunik).

#Example of a sharp RDD
#let's prepare example data
x <- runif(1000, -2, 2)
y <- 3 + 2 * x + 5 * (x>=1) + rnorm(1000) #note the condition

plot(x,y,
     col = "blue",
     cex = 0.35)
lines(c(1, 1), c(-20, 30), col = "black", lty = 2, lwd = 2)

#step 1: create an rdd_data object (x = running variable, cutpoint = threshold value)
dane_srdd <- rdd_data(y = y, x = x, cutpoint = 1) #VERY IMPORTANT - define the cutpoint!
summary(dane_srdd) #note the "Type: Sharp"
plot(dane_srdd)    #binned plot

#step 2a: non-parametric estimation (local linear regression)
bw <- rdd_bw_ik(dane_srdd) #optimal Imbens-Kalyanaraman bandwidth
bw
srdd_np <- rdd_reg_np(dane_srdd, bw = bw)
summary(srdd_np) #LATE = local average treatment effect
plot(srdd_np)

#step 2b: parametric estimation (a polynomial of a given order)
srdd_lm <- rdd_reg_lm(dane_srdd, order = 1)
summary(srdd_lm)
plot(srdd_lm)

#McCrary test - whether observations do not "bunch up" just past the threshold
#(if units could manipulate the running variable, RDD would be invalid)
dens_test(dane_srdd)


#Example of a fuzzy RDD
set.seed(2020)
x <- runif(1000, -2, 2)
S<-rbinom(1000, 1, prob = 0.8)
y <- 3 + 2 * x + 5 *S* (x>=1) + rnorm(1000,0,2) #note the condition

plot(x,y,
     col = "blue",
     cex = 0.35)
lines(c(1, 1), c(-20, 30), col = "black", lty = 2, lwd = 2)

#S = actual treatment status: no one below the threshold, 80% of observations above it
S[x < 1] <- 0

#in a fuzzy RDD we additionally supply z = a variable describing actual treatment
dane_frdd <- rdd_data(y = y, x = x, z = S, cutpoint = 1)
summary(dane_frdd) #now "Type: Fuzzy"

frdd <- rdd_reg_np(dane_frdd)
summary(frdd) #LATE = local average treatment effect

#for comparison: let's treat the same data as sharp RDD
dane_srdd2 <- rdd_data(y = y, x = x, cutpoint = 1)
srdd2 <- rdd_reg_np(dane_srdd2)
summary(srdd2)

res<-c(fuzzy = rdd_coef(frdd), sharp = rdd_coef(srdd2))
res #the sharp-variant effect is understated - since only 80% of obs. above the threshold were actually treated


#####################################
#4. The synthetic control method - the BLIK example
#Data: 22 countries (Poland + 21 donor pool countries), 2000-2024
#Outcome variable: the value of electronic payments as a % of household consumption expenditure

dane <- read.csv("W3_BLIK_SCM.csv")
str(dane)
table(dane$Country)

#it's always worth looking at the data before estimation
plot(ELECTRONIC_2_CONS ~ Year, data = dane[dane$Country == "Poland", ],
     type = "l", col = "blue", lwd = 2, ylim = c(0, 60),
     ylab = "Electronic payments / HFCE (%)", xlab = "Year")
for (k in setdiff(unique(dane$Country), "Poland"))
  lines(ELECTRONIC_2_CONS ~ Year, data = dane[dane$Country == k, ], col = "grey80")
lines(ELECTRONIC_2_CONS ~ Year, data = dane[dane$Country == "Poland", ],
      col = "blue", lwd = 2)
abline(v = 2015, lty = 2)

#the whole process in a single pipe
blik <- dane |>
  synthetic_control(outcome = ELECTRONIC_2_CONS,
                    unit = Country,
                    time = Year,
                    i_unit = "Poland",   #the treated unit
                    i_time = 2015,       #the moment of intervention
                    generate_placebos = TRUE) |>   #generate placebos right away, for inference

  #predictors: averages over the PRE-treatment period
  generate_predictor(time_window = 2000:2014,
                     pkb_pc      = mean(GDP_PER_CAPITA_TH, na.rm = TRUE),
                     rozwoj_fin  = mean(FIN_DEV, na.rm = TRUE),
                     internet    = mean(INTERNET_ACCESS, na.rm = TRUE),
                     edukacja    = mean(HIGHER_EDUCATION, na.rm = TRUE),
                     rzady_prawa = mean(RULE_OF_LAW, na.rm = TRUE)) |>

  #weight selection: minimising RMSPE in the pre-treatment period
  generate_weights(optimization_window = 2000:2014) |>

  generate_control()

#1) pre-treatment fit and the post-treatment gap
blik |> plot_trends()
blik |> plot_differences()

#2) country weights and predictor weights
blik |> plot_weights()
blik |> grab_unit_weights() |> subset(weight > 0.001)

#3) predictor balance: Poland vs. Synthetic Poland
blik |> grab_balance_table()

#4) inference: placebo tests and the permutation p-value
blik |> plot_placebos()                    #every country treated as "falsely treated"
blik |> plot_placebos(prune = TRUE)        #excluding countries with a poor pre-treatment fit
blik |> plot_mspe_ratio()                  #post/pre RMSPE ratio
blik |> grab_significance()                #permutation p-value

#5) the gap itself, year by year
luka <- blik |> grab_synthetic_control()
luka$gap <- luka$real_y - luka$synth_y
luka

#Published results (Stata, allsynth) for comparison:
#weights: Hungary 0.388, Latvia 0.359, Bulgaria 0.124, Albania 0.071, Greece 0.058
#pre-treatment RMSPE: 0.24 pp; 2023 gap: 5.8 pp; 2024 gap: 8.9 pp
#NOTE: optimisation of the V matrix is often weakly identified, so R may return
#somewhat different weights than Stata. Check whether the qualitative conclusions still hold.


#####################################
# Homework
#1. For the fuzzy RDD data, check how the LATE estimate changes when you change the bandwidth (bw).
#2. In the staggered DiD example, change the effect strength so that it is IDENTICAL across all
#   cohorts. Does TWFE then recover the true effect?
#3. Repeat the SCM estimation while removing Hungary from the donor pool (leave-one-out).
#   How does the estimated gap change?
#4. Repeat the SCM estimation using 2012 as the moment of intervention (a placebo in time).
#   Does the gap appear before 2015?