Fix Y permutation error, implement placebo SEs

Author: nilshg

src/SyntheticDiD.jl

@@ -12,6 +12,10 @@ a unit weight and λ is a time period weight.
 
 The implementation follows the author's reference implementation in the R package `synthdid`.
 
+NOTE: The implementation assumes that the outcomes and treatment matrices are sorted such that
+treated units come last in both `Y` and `W`. It therefore checks whether `W` is sorted and swaps
+rows in both `Y` and `W` to sort both matrices accordingly if required. 
+
 """
 struct SyntheticDiD{T1}# <: SCM where T1
     tp::T1
@@ -35,27 +39,63 @@ function isfitted(x::SyntheticDiD)
     all(!iszero, x.τ̂) || all(!iszero, x.se_τ̂)
 end
 
-function fit!(x::SyntheticDiD{BalancedPanel{SingleUnitTreatment{Continuous}}})
+function fit!(x::SyntheticDiD{BalancedPanel{SingleUnitTreatment{Continuous}}};
+    se = nothing)
     (;Y, W) = x.tp
 
     T₀ = length_T₀(x.tp)
     N₁ = 1 # length(treated_ids(x)) - has to be one given type of BalancedPanel in this method
     N₀ = size(Y, 1) - N₁
 
     # Implementation of synthdid below assumes that treated unit comes last
-    swaprows!(Y, treated_ids(x.tp), size(Y, 1))
+    if !issorted(any(x) for x ∈ eachrow(W)) 
+        treated_row = treated_ids(x.tp)
+        swaprows!(Y, treated_row, size(Y, 1))
+        swaprows!(W, treated_row, size(W, 1))
+    end
 
     # Get estimate and store in SyntheticDiD object
     estimate = synthdid_estimate(Y, N₀, T₀)   
     x.τ̂ .= estimate.τ̂
     x.λ̂ .= estimate.λ̂
     x.ω̂ .= estimate.ω̂
 
+    if !isnothing(se)
+        if se != :placebo
+            error("Only placebo standard error estimation is currently implemented")
+        else
+            V̂ = get_V̂_τ(x).V̂
+            x.se_τ̂ .= √(V̂)
+        end
+    end
+
     # return SynhtDiD object
     return x
 end
 
 
+function get_V̂_τ(x::SyntheticDiD{BalancedPanel{SingleUnitTreatment{Continuous}}})
+    (;Y, W) = x.tp
+
+    Y_untreated = Y[Not(treated_ids(x.tp)), :]
+
+    T₀ = length_T₀(x.tp)
+    N₁ = 1
+    N₀ = size(Y_untreated, 1) - N₁
+
+    τ̂s = zeros(N₀)
+
+    for i ∈ 1:N₀-1
+        swaprows!(Y_untreated, i, size(Y_untreated, 1))
+        τ̂s[i] = synthdid_estimate(Y_untreated, N₀, T₀).τ̂
+    end
+
+    swaprows!(Y_untreated, 1, size(Y_untreated, 1))
+    τ̂s[end] = synthdid_estimate(Y_untreated, N₀, T₀).τ̂
+
+    return (V̂ = var(τ̂s), τ̂s = τ̂s)
+end
+
 function synthdid_estimate(Y, N₀, T₀;
     ω_intercept = true, λ_intercept = true,
     max_iter = 10_000, max_iter_pre_sparsify = 100, sparsify_function = sparsify_weights!)
@@ -196,7 +236,7 @@ end
 
 
 function swaprows!(X::AbstractMatrix, i::Integer, j::Integer)
-    @inbounds for k = 1:size(X,2)
+    @inbounds for k = axes(X,2)
         X[i, k], X[j, k] = X[j, k], X[i, k]
     end
 end

test/runtests.jl

@@ -39,6 +39,12 @@ end
   fit!(sdid_model)
 
   @test only(sdid_model.τ̂) ≈ -15.604 atol=0.1
+
+  @test_throws "placebo" fit!(sdid_model, se = :jackknife)
+
+  fit!(sdid_model, se = :placebo)
+
+  @test only(sdid_model.se_τ̂) ≈ 9.3 atol=0.1
 end