implement dd estimator

Author: jbytecode

CHANGELOG.md

@@ -1,5 +1,7 @@
 # v0.1.7 (Upcoming Release)
 
+- Implement IV (Instrumental Variables) regression estimator
+- Implement DD (Dagenais and Dagenais, 1997) moments estimator
 
 # v0.1.6
 

Project.toml

@@ -1,6 +1,6 @@
 name = "ErrorsInVariables"
 uuid = "faeb8267-8352-4c3e-81c7-04150084d5bf"
-version = "0.1.6"
+version = "0.1.7"
 authors = ["jbytecode <mhsatman@gmail.com>"]
 
 [deps]

src/dd.jl

@@ -17,13 +17,17 @@ Computes the Higher Moment Estimator for a linear regression model with errors i
 - Dagenais, Marcel G., and Denyse L. Dagenais. "Higher moment estimators for linear regression models with errors in the variables." 
 Journal of Econometrics 76.1-2 (1997): 193-221.
 """
-function dd(xvec::Vector, yvec::Vector)::Matrix 
+function dd(xmat::Matrix, yvec::Vector)::Matrix 
 
-    n = length(xvec)
+    n, p = size(xmat)
 
-    IK = ones(n)
+    o = ones(p)
+    IK = I(p)
 
-    x = xvec .- mean(xvec)
+    x = similar(xmat)
+    for i in 1:p
+        x[:, i] = xmat[:, i] .- mean(xmat[:, i])
+    end
 
     y = yvec .- mean(yvec)
 
@@ -33,11 +37,15 @@ function dd(xvec::Vector, yvec::Vector)::Matrix
 
     z3 = y .* y
 
-    z4 = x .* x .* x .- 3 .* x .* (mean(x'x / n) .* IK)
+    z4 = x .* x .* x .- 3 .* x * (mean(x'x ./ n) * IK)
 
-    # z5 = x .* x .* y .- 2 .* x .* (mean(y'x / n) .* IK) .- y * (IK' * mean(x'x / n) .* IK)
+    z5 = x .* x .* y .- 2 .* x * (mean(y'x ./ n) .* IK) .- y .* (mean(x'x / n) .* ones(n))
 
-    return hcat(z1, z2, z3, z4)
+    z6 = x .* y .* y .- x * (mean(y'y ./ n) .* IK) .- 2 * y .* (mean(y'x / n) .* ones(n))
+
+    z7 = y .* y .* y .- 3 .* y .* (mean(y'y ./ n))
+
+    return hcat(z1, z2, z3, z4, z5, z6, z7)
 
 end 
 

src/iv.jl

@@ -18,7 +18,7 @@ computed using the formula:
 If the number of instruments is greater than the number of predictors, the IV estimator is
 computed using the formula:
 ```julia
-    part = X' * Z * inv(Z' * Z) * Z' * X
+    part = X' * Z * inv(Z' * Z) * Z' 
     betas = inv(part * X) * part * y
 ```
 """
@@ -33,7 +33,7 @@ function iv(X::Matrix, y::Vector, Z::Matrix)::Vector
     if numberofinstruments == numberofpredictors
         return inv(Z'X)Z'y
     else
-        part = X' * Z * inv(Z' * Z) * Z' * X
+        part = X' * Z * inv(Z' * Z) * Z' 
         betas = inv(part * X) * part * y
         return betas
     end

test/runtests.jl

@@ -7,12 +7,12 @@ function euclidean(v1::Vector, v2::Vector)
 end
 
 include("testdd.jl")
-#include("testiv.jl")
-#include("testcga.jl")
-#include("testestimator.jl")
-#include("testmeive.jl")
-#include("testorthogonalregression.jl")
-#include("testdeming.jl")
-#include("testsimex.jl")
+include("testiv.jl")
+include("testcga.jl")
+include("testestimator.jl")
+include("testmeive.jl")
+include("testorthogonalregression.jl")
+include("testdeming.jl")
+include("testsimex.jl")
 
 

test/testdd.jl

@@ -1,18 +1,46 @@
 using Test 
 using ErrorsInVariables
+using Random
 
 
 @testset "DD estimator" verbose=true begin 
 
     @testset "Simple Example" verbose=true begin 
 
-        X = [1.0, 2.0, 3.0, 4.0, 5.0, 6, 7, 8]
+        eps = 0.0001
 
-        y = [2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0]
+        rng = Random.MersenneTwister(1234)
 
-        result = dd(X, y)
+        n = 30
+
+        deltax = randn(rng, n) * sqrt(3.0)
+
+        cleanx = randn(rng, n) * sqrt(7.0)
+        cleanx2 = randn(rng, n) * sqrt(7.0)
+        cleanx3 = randn(rng, n) * sqrt(7.0)
+
+        e = randn(rng, n) * sqrt(5.0)
+
+        y = 20.0 .+ 10.0 .* cleanx .+ 15.0 .* cleanx2 .+ 13.0 .* cleanx3 .+ e
+
+        dirtyx = cleanx .+ deltax
+
+        Xd = hcat(ones(n), dirtyx, cleanx2, cleanx3)
+        Xc = hcat(ones(n), cleanx, cleanx2, cleanx3)
+
+        dirtybetas = Xd \ y
+
+        cleanbetas = Xc \ y
+
+        Z = dd(hcat(dirtyx, cleanx2, cleanx3), y)
+
+        result = iv(Xd, y, Z)
+
+        expected_dd_betas = [18.488555621677015, 6.283197918953448, 15.765140901269259, 12.061208813549317]
+
+        for i in eachindex(expected_dd_betas)
+            @test isapprox(result[i], expected_dd_betas[i], atol=eps)
+        end
 
-        println("-----------------------")
-        display(result)
     end 
 end