rm(list=ls()) #czyszczenie obszaru roboczego set.seed(2020) #aby wszyscy otrzymali takie same wyniki #Ustaw katalog roboczy na folder z plikiem W3_BLIK_SCM.csv, np.: #W RStudio: Session -> Set Working Directory -> To Source File Location #uwaga: skrypt korzysta z operatora potoku |> (wymaga R w wersji 4.1 lub nowszej) #uwaga: require() samo nie ładuje pakietu po instalacji - stąd library() w nawiasie klamrowym if(!require("rddtools")) {install.packages("rddtools"); library("rddtools")} if(!require("tidysynth")) {install.packages("tidysynth"); library("tidysynth")} ##################################### #1. Difference-in-differences - ilustracja #tworzymy wykres dla grupy kontrolnej plot(c(0, 1), c(5, 7), type = "p", ylim = c(4, 12), xlim = c(-0.2, 1.2), main = "Estymator DID", xlab = "Czas", ylab = "y", col = "grey", pch = 19, cex = 1.2, xaxt = "n", yaxt = "n") axis(1, at = c(0, 1), labels = c("przed", "po")) #grupa poddana oddziaływaniu points(c(0, 1), c(7, 11), col = "blue", pch = 19, cex = 1.2) #grupa poddana oddziaływaniu gdyby zachowała ten sam trend co grupa kontrolna points(1,9, col = "lightblue", pch = 19, cex = 1.2) #łączymy punkty 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) #dodajemy opis text(1, 9.7, expression(hat(beta)[1]^{DID}), cex = 0.8, pos = 4) text(1, 10.3, "efekt oddziaływania", cex = 0.8, pos = 4) #a teraz przykład na wygenerowanych danych N<-250 #ustalamy wielkość próby #określamy efekt oddziaływania TreatEff<-5 #tworzymy zmienną która podzieli N obserwacji na grupy S<- c(rep(0, N/2), rep(1, N/2)) #tworzymy dane przed i po oddziaływaniu #uwaga na nawias: 1:(N/2) to obserwacje 1-125, a 1:N/2 to (1:N)/2 - czyli zupełnie co innego! 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 #policzmy efekt mean(y_po[S == 1]) - mean(y_przed[S == 1]) - (mean(y_po[S == 0]) - mean(y_przed[S == 0])) #dlaczego nie wyszło nam dokładnie 5? #policzmy za pomocą KMNK summary(lm(y_po - y_przed ~ S)) ##################################### #2. DiD przy zróżnicowanych momentach wejścia w życie (staggered adoption) #Pokazujemy, dlaczego zwykła regresja TWFE może zawieść, gdy interwencja #wchodzi w życie w różnych momentach, a efekty są zróżnicowane. N_j <- 300; T_max <- 20 panel <- expand.grid(id = 1:N_j, t = 1:T_max) #trzy kohorty wchodzą w życie w różnych momentach kohorta <- sample(c(5, 10, 15), N_j, replace = TRUE) panel$g <- kohorta[panel$id] #efekt PRAWDZIWY: rośnie z czasem od momentu wejścia i różni się między kohortami #(im wcześniejsza kohorta, tym silniejszy efekt - to wystarczy, by zaburzyć 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)) #prawdziwy średni efekt wśród obserwacji objętych działaniem prawdziwy <- mean(panel$efekt[panel$D == 1]) #estymator TWFE (efekty stałe jednostki i czasu) twfe <- lm(y ~ D + factor(id) + factor(t), data = panel) oszacowany <- coef(twfe)["D"] c(prawdziwy = prawdziwy, TWFE = oszacowany) #TWFE nie odtwarza prawdziwego efektu - część porównań używa jednostek #JUŻ objętych działaniem jako grupy kontrolnej (Goodman-Bacon 2021) #Rozwiązania: pakiety did (Callaway, Sant'Anna), fixest::sunab (Sun, Abraham), #didimputation. Zob. też wykres event study: fixest::feols + iplot ##################################### #3. Regression discontinuity design #przykład: o przyjęciu na studia decyduje wynik testu (próg 70%) x <- runif(1000, 0, 1) y <- as.numeric(x >= 0.7) #odpowiednik pętli: 1 gdy próg przekroczony, 0 w przeciwnym razie plot(x,y, col = "blue", cex = 0.35, xlab = "Wynik testu", ylab = "Przyjęcie na studia") #Uwaga techniczna: pakiet rdd (funkcja RDestimate) został zdjęty z CRAN 10.07.2025. #Korzystamy z pakietu rddtools, który wrócił na CRAN 29.10.2025 (wersja 2.0.2) #i oferuje ZARÓWNO estymację nieparametryczną, JAK I parametryczną. #Alternatywa stosowana dziś w literaturze: pakiet rdrobust (Calonico, Cattaneo, Titiunik). #Przykład ostrego RDD (sharp RDD) #Przygotujemy przykładowe dane x <- runif(1000, -2, 2) y <- 3 + 2 * x + 5 * (x>=1) + rnorm(1000) #zwróćcie uwagę na warunek plot(x,y, col = "blue", cex = 0.35) lines(c(1, 1), c(-20, 30), col = "black", lty = 2, lwd = 2) #krok 1: tworzymy obiekt rdd_data (x = zmienna progowa, cutpoint = wartość progu) dane_srdd <- rdd_data(y = y, x = x, cutpoint = 1) #BARDZO WAŻNE - zdefiniować cutpoint! summary(dane_srdd) #zwróćcie uwagę na "Type: Sharp" plot(dane_srdd) #wykres z uśrednionymi przedziałami (binned plot) #krok 2a: estymacja nieparametryczna (lokalna regresja liniowa) bw <- rdd_bw_ik(dane_srdd) #optymalne pasmo Imbensa-Kalyaramana bw srdd_np <- rdd_reg_np(dane_srdd, bw = bw) summary(srdd_np) #LATE = local average treatment effect plot(srdd_np) #krok 2b: estymacja parametryczna (wielomian zadanego rzędu) srdd_lm <- rdd_reg_lm(dane_srdd, order = 1) summary(srdd_lm) plot(srdd_lm) #test McCrary'ego - czy obserwacje nie "gromadzą się" tuż za progiem #(gdyby jednostki mogły manipulować zmienną progową, RDD byłoby nieważne) dens_test(dane_srdd) #Przykład nieostrego RDD (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) #zwróćcie uwagę na warunek plot(x,y, col = "blue", cex = 0.35) lines(c(1, 1), c(-20, 30), col = "black", lty = 2, lwd = 2) #S = faktyczne poddanie oddziaływaniu: poniżej progu nikt, powyżej progu 80% obserwacji S[x < 1] <- 0 #w fuzzy RDD podajemy dodatkowo z = zmienną opisującą faktyczne oddziaływanie dane_frdd <- rdd_data(y = y, x = x, z = S, cutpoint = 1) summary(dane_frdd) #teraz "Type: Fuzzy" frdd <- rdd_reg_np(dane_frdd) summary(frdd) #LATE = local average treatment effect #dla porównania: potraktujmy te same dane jak 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 #efekt w wariancie sharp jest zaniżony - bo tylko 80% obs. powyżej progu było objętych działaniem ##################################### #4. Metoda kontroli syntetycznej - przykład BLIK #Dane: 22 kraje (Polska + 21 krajów donor pool), lata 2000-2024 #Zmienna objaśniana: wartość płatności elektronicznych jako % spożycia gosp. domowych dane <- read.csv("W3_BLIK_SCM.csv") str(dane) table(dane$Country) #zawsze warto zobaczyć dane przed estymacją plot(ELECTRONIC_2_CONS ~ Year, data = dane[dane$Country == "Poland", ], type = "l", col = "blue", lwd = 2, ylim = c(0, 60), ylab = "Płatności elektroniczne / HFCE (%)", xlab = "Rok") 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) #cały proces w jednym potoku (pipe) blik <- dane |> synthetic_control(outcome = ELECTRONIC_2_CONS, unit = Country, time = Year, i_unit = "Poland", #jednostka objęta interwencją i_time = 2015, #moment interwencji generate_placebos = TRUE) |> #od razu licz placebo dla wnioskowania #predyktory: średnie z okresu PRZED interwencją 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)) |> #dobór wag: minimalizacja RMSPE w okresie przed interwencją generate_weights(optimization_window = 2000:2014) |> generate_control() #1) dopasowanie przed interwencją i luka po interwencji blik |> plot_trends() blik |> plot_differences() #2) wagi krajów i wagi predyktorów blik |> plot_weights() blik |> grab_unit_weights() |> subset(weight > 0.001) #3) bilans predyktorów: Polska vs Syntetyczna Polska blik |> grab_balance_table() #4) wnioskowanie: testy placebo i p-value permutacyjne blik |> plot_placebos() #wszystkie kraje jako "fałszywie objęte" blik |> plot_placebos(prune = TRUE) #bez krajów o słabym dopasowaniu przed interwencją blik |> plot_mspe_ratio() #iloraz RMSPE po/przed interwencją blik |> grab_significance() #p-value permutacyjne #5) sama luka rok po roku luka <- blik |> grab_synthetic_control() luka$gap <- luka$real_y - luka$synth_y luka #Wyniki opublikowane (Stata, allsynth) dla porównania: #wagi: Węgry 0.388, Łotwa 0.359, Bułgaria 0.124, Albania 0.071, Grecja 0.058 #RMSPE przed interwencją: 0.24 pp; luka 2023: 5.8 pp; luka 2024: 8.9 pp #UWAGA: optymalizacja macierzy V bywa słabo zidentyfikowana, więc R może zwrócić #nieco inne wagi niż Stata. Sprawdźcie, czy wnioski jakościowe pozostają takie same. ##################################### # Praca domowa #1. Dla danych fuzzy RDD sprawdź, jak zmienia się oszacowanie LATE, gdy zmienisz pasmo (bw). #2. W przykładzie staggered DiD zmień siłę efektu tak, aby była JEDNAKOWA we wszystkich # kohortach. Czy TWFE odtwarza wtedy prawdziwy efekt? #3. Powtórz estymację SCM, usuwając Węgry z donor pool (leave-one-out). # Jak zmienia się oszacowana luka? #4. Powtórz estymację SCM, przyjmując jako moment interwencji rok 2012 (placebo w czasie). # Czy luka pojawia się przed 2015 rokiem?