# Chatpter 5 (cont.). Advanced Causal Inference in Julia ```julia import Pkg; Pkg.add(["DiffinDiffs", "RegressionDiscontinuity", "SynthControl", "Dates", "Plots"]); ``` ```  Resolving package versions...  No Changes to `~/.julia/environments/v1.11/Project.toml`  No Changes to `~/.julia/environments/v1.11/Manifest.toml` ``` ## Instrumental Variable ```julia using DataFrames, RDatasets, FixedEffectModels df = dataset("plm", "Cigar") first(df, 5) ``` 5×9 DataFrame | Row | State | Year | Price | Pop | Pop16 | CPI | NDI | Sales | Pimin | | --- | ----- | ----- | ------- | ------- | ------- | ------- | ------- | ------- | ------- | | | Int64 | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | | 1 | 1 | 63 | 28.6 | 3383.0 | 2236.5 | 30.6 | 1558.3 | 93.9 | 26.1 | | 2 | 1 | 64 | 29.8 | 3431.0 | 2276.7 | 31.0 | 1684.07 | 95.4 | 27.5 | | 3 | 1 | 65 | 29.8 | 3486.0 | 2327.5 | 31.5 | 1809.84 | 98.5 | 28.9 | | 4 | 1 | 66 | 31.5 | 3524.0 | 2369.7 | 32.4 | 1915.16 | 96.4 | 29.5 | | 5 | 1 | 67 | 31.6 | 3533.0 | 2393.7 | 33.4 | 2023.55 | 95.5 | 29.6 | For example, if we have an instrumental variable `Price` predicting `CPI` which predicts `Sales` (theoretically incorrect, just for illustration), the formula here should be: ```julia # DV ~ exogenous variables + (endogenous variables ~ instrumental variables) + fe(fixedeffect variable) reg(df, @formula(Sales ~ NDI + (CPI ~ Price ) + fe(State))) ``` ``` FixedEffectModel ======================================================================= Number of obs: 1380 Converged: true dof (model): 2 dof (residuals): 1333 R²: -0.327 R² adjusted: -0.373 F-statistic: 42.695 P-value: 0.000 F-statistic (first stage): 39.4447 P-value (first stage): 0.000 R² within: -3.259 Iterations: 1 ======================================================================= Estimate Std. Error t-stat Pr(>|t|) Lower 95% Upper 95% ─────────────────────────────────────────────────────────────────────── NDI 0.0271309 0.00634835 4.2737 <1e-04 0.0146771 0.0395848 CPI -3.71621 0.817692 -4.54475 <1e-05 -5.32031 -2.1121 ======================================================================= ``` ## Difference-in-difference ```julia using DiffinDiffs hrs = DiffinDiffsBase.exampledata("hrs") ``` ``` 3280×11 VecColumnTable:  Row │ hhidpn  wave  wave_hosp  oop_spend  riearnsemp  rwthh  male  spouse  ⋯ │ Int64  Int64  Int64  Float64  Float64  Int64  Int64  Int64  ⋯ ──────┼───────────────────────────────────────────────────────────────────────── 1 │ 1 10 10 6532.91 6.37159e5 4042 0 0 ⋯ 2 │ 1 8 10 1326.93 3.67451e5 3975 0 0 ⋯ 3 │ 1 11 10 1050.33 74130.5 3976 0 0 ⋯ 4 │ 1 9 10 979.418 84757.4 3703 0 0 ⋯ 5 │ 1 7 10 5498.68 1.66128e5 5295 0 0 ⋯ 6 │ 2 8 8 41504.0 0.0 5187 0 1 ⋯ 7 │ 2 7 8 3672.86 0.0 4186 0 1 ⋯ 8 │ 2 10 8 1174.19 0.0 3729 0 1 ⋯ 9 │ 2 11 8 6909.59 0.0 3906 0 1 ⋯ 10 │ 2 9 8 1130.1 0.0 5453 0 1 ⋯ 11 │ 3 9 10 2773.45 8475.74 4419 0 1 ⋯ ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 3270 │ 654 10 10 2876.32 45511.3 4372 0 0 ⋯ 3271 │ 655 7 9 1910.74 25476.5 7682 0 1 ⋯ 3272 │ 655 11 9 12402.5 34721.6 9203 0 1 ⋯ 3273 │ 655 8 9 1530.0 45000.0 8461 0 1 ⋯ 3274 │ 655 9 9 7373.89 10359.2 9345 0 1 ⋯ 3275 │ 655 10 9 673.568 38229.5 8420 0 1 ⋯ 3276 │ 656 11 8 3020.78 0.0 1930 0 0 ⋯ 3277 │ 656 8 8 2632.0 0.0 4810 0 0 ⋯ 3278 │ 656 9 8 657.34 0.0 4768 0 0 ⋯ 3279 │ 656 10 8 782.795 0.0 1909 0 0 ⋯ 3280 │ 656 7 8 4182.39 0.0 4374 0 0 ⋯ ``` ```julia r = @did(Reg, data = hrs, dynamic(:wave, -1), notyettreated(11), vce = Vcov.cluster(:hhidpn), yterm = term(:oop_spend), treatname = :wave_hosp, treatintterms = (), xterms = (fe(:wave) + fe(:hhidpn))) ``` ``` ────────────────────────────────────────────────────────────────────── Summary of results: Regression-based DID ────────────────────────────────────────────────────────────────────── Number of obs: 2624 Degrees of freedom: 14 F-statistic: 6.42 p-value: <1e-07 ────────────────────────────────────────────────────────────────────── Cohort-interacted sharp dynamic specification ────────────────────────────────────────────────────────────────────── Number of cohorts: 3 Interactions within cohorts: 0 Relative time periods: 5 Excluded periods: -1 ────────────────────────────────────────────────────────────────────── Fixed effects: fe_hhidpn fe_wave ────────────────────────────────────────────────────────────────────── Converged: true Singletons dropped: 0 ────────────────────────────────────────────────────────────────────── ``` ```julia a = agg(r, :rel) ``` ``` ─────────────────────────────────────────────────────────────────── Estimate Std. Error t Pr(>|t|) Lower 95% Upper 95% ─────────────────────────────────────────────────────────────────── rel: -3 591.046 1273.08 0.46 0.6425 -1905.3 3087.39 rel: -2 352.639 697.78 0.51 0.6133 -1015.62 1720.9 rel: 0 2960.04 540.989 5.47 <1e-07 1899.23 4020.86 rel: 1 529.767 586.831 0.90 0.3667 -620.935 1680.47 rel: 2 800.106 1010.81 0.79 0.4287 -1181.97 2782.18 ─────────────────────────────────────────────────────────────────── ``` ## Regression Discontinuity ```julia using RegressionDiscontinuity data = RDData(RegressionDiscontinuity.Lee08()) ``` ``` RDData{Vector{Float64}, RunningVariable{Float64, Float64, Vector{Float64}}}([0.0, 0.0, 0.0, 0.251144230365753, 0.23875193297862998, 0.267181247472763, 0.262862056493759, 0.0, 0.34071296453476, 0.6447104215621949 … 0.467111736536026, 0.712594926357269, 0.718958139419556, 0.791856050491333, 0.8541748523712159, 0.677568554878235, 0.808446884155273, 0.505452692508698, 1.0, 0.760513544082642], [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) ``` ```julia fit(NaiveLocalLinearRD(kernel=Rectangular(), bandwidth=ImbensKalyanaraman()), data.ZsR, data.Ys) ``` ``` Local linear regression for regression discontinuity design ⋅⋅⋅⋅ Naive inference (not accounting for bias) ⋅⋅⋅⋅ Rectangular kernel (U[-0.5,0.5]) ⋅⋅⋅⋅ Imbens Kalyanaraman bandwidth ⋅⋅⋅⋅ Eicker White Huber variance ──────────────────────────────────────────────────────────────────────────────────────────────── h τ̂ se bias z p-val Lower 95% Upper 95% ──────────────────────────────────────────────────────────────────────────────────────────────── Sharp RD estimand 0.462024 0.08077 0.0087317 unaccounted 9.25 <1e-99 0.0636562 0.0978838 ──────────────────────────────────────────────────────────────────────────────────────────────── ``` ## Synthetic Control ```julia using SynthControl, Dates df = load_brexit() ``` 897×3 DataFrame872 rows omitted | Row | country | quarter | realgdp | | --- | -------------- | ---------- | ------- | | | String | Date | Float64 | | 1 | Australia | 2009-01-01 | 1.04 | | 2 | Austria | 2009-01-01 | -1.53 | | 3 | Belgium | 2009-01-01 | -1.15 | | 4 | Canada | 2009-01-01 | -2.28 | | 5 | Denmark | 2009-01-01 | -1.42 | | 6 | Finland | 2009-01-01 | -6.8 | | 7 | France | 2009-01-01 | -1.67 | | 8 | Germany | 2009-01-01 | -4.49 | | 9 | Greece | 2009-01-01 | -4.74 | | 10 | Iceland | 2009-01-01 | -7.61 | | 11 | Ireland | 2009-01-01 | -0.34 | | 12 | Italy | 2009-01-01 | -2.75 | | 13 | Japan | 2009-01-01 | -4.88 | | ⋮ | ⋮ | ⋮ | ⋮ | | 886 | Italy | 2018-07-01 | -0.82 | | 887 | Japan | 2018-07-01 | 9.54 | | 888 | Luxembourg | 2018-07-01 | 28.41 | | 889 | Netherlands | 2018-07-01 | 10.69 | | 890 | New Zealand | 2018-07-01 | 24.06 | | 891 | Norway | 2018-07-01 | 12.82 | | 892 | Portugal | 2018-07-01 | 2.23 | | 893 | Spain | 2018-07-01 | 5.74 | | 894 | Sweden | 2018-07-01 | 22.48 | | 895 | Switzerland | 2018-07-01 | 14.35 | | 896 | United Kingdom | 2018-07-01 | 15.72 | | 897 | United States | 2018-07-01 | 19.32 | ```julia bp = BalancedPanel(df, "United Kingdom" => Date(2016, 7, 1); id_var = :country, t_var = :quarter, outcome_var = :realgdp) ``` ``` Balanced Panel - single treated unit, continuous treatment Treated unit: United Kingdom Number of untreated units: 22 First treatment period: 2016-07-01 Number of pretreatment periods: 30 Number of treatment periods: 9 ``` ```julia s = SimpleSCM(bp) ``` ``` Synthetic Control Model Treatment panel: Model is not fitted Balanced Panel - single treated unit, continuous treatment Treated unit: United Kingdom Number of untreated units: 22 First treatment period: 2016-07-01 Number of pretreatment periods: 30 Number of treatment periods: 9 ``` ```julia fit!(s) ``` ``` Synthetic Control Model Treatment panel: Model is fitted Impact estimates: [-0.54, -0.31, -0.206, -0.732, -1.241, -1.482, -1.818, -2.327, -1.994] Balanced Panel - single treated unit, continuous treatment Treated unit: United Kingdom Number of untreated units: 22 First treatment period: 2016-07-01 Number of pretreatment periods: 30 Number of treatment periods: 9 ``` ```julia using Plots plot(s) ``` ![](/files/3msIh3v75hyhQ4Uv61Bk) ## SyntheticDiD ```julia sp = load_smoking_panel() ``` ``` Balanced Panel - single treated unit, continuous treatment Treated unit: 3 Number of untreated units: 38 First treatment period: 1989 Number of pretreatment periods: 19 Number of treatment periods: 12 ``` ```julia sdid_model = SyntheticDiD(sp) ``` ``` Synthetic Difference-in-Differences Model Model is not fitted ``` ```julia fit!(sdid_model) ``` ``` Synthetic Difference-in-Differences Model Model is fitted Impact estimate: -15.604 ``` --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://data-julia.rongxin.me/5.2.models.jl.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.