gragusa / InstrumentalVariables.jl

A Julia package for Instrumental Variables models
Other
3 stars 1 forks source link

add standard TSLS standard errors #2

Open remlapmot opened 8 years ago

remlapmot commented 8 years ago

It's great that you have interfaced to the heteroskedasticity robust standard errors but would you mind adding the standard TSLS SEs. For the example you're using from the AER R package the default SEs from the ivreg() output are as follows:

library(Ecdat)
library(AER)
data(Cigarette)
attach(Cigarette)
y <- log(packpc)
x1 <- log(income/pop/cpi)
x2 <- log(avgprs/cpi)
z1 <- (taxs - tax)/cpi
z2 <- tax/cpi
iv <- ivreg(y ~ x1 + x2 | x1 + z1 + z2, subset=year==1995)
summary(iv)
Call:
ivreg(formula = y ~ x1 + x2 | x1 + z1 + z2, subset = year == 
    1995)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.6006931 -0.0862222 -0.0009999  0.1164699  0.3734227 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.8950     1.0586   9.348 4.12e-12 ***
x1            0.2804     0.2386   1.175    0.246    
x2           -1.2774     0.2632  -4.853 1.50e-05 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1879 on 45 degrees of freedom
Multiple R-Squared: 0.4294, Adjusted R-squared: 0.4041 
Wald test: 13.28 on 2 and 45 DF,  p-value: 2.931e-05 

cheers, Tom

gragusa commented 8 years ago

Thank you for pointing this out. Code for the homoskedastic standard errors was written by there was no way to call it. I have changed things so that the default show method gives these standard errors as default

using RDatasets
using DataFramesMeta
cig = dataset("Ecdat", "Cigarette")
cig[:lrincome] = @with(cig, log(:Income./:POP./:CPI));
cig[:lrprice] = @with(cig, log(:AvgPrs./:CPI));
cig[:tdiff] = @with(cig, (:Taxs - :Tax)./:CPI)
cig[:rtax] = @with(cig, :Tax./:CPI)
cig95 = @ix(cig, :Year .== 1995)
y = convert(Array, log(cig95[:PackPC]));
x = convert(Array, [ones(size(cig95,1)) cig95[:lrincome] cig95[:lrprice]])
z = convert(Array, [ones(size(cig95,1)) cig95[:lrincome] cig95[:tdiff] cig95[:rtax]])

iv = ivreg(x, z, y)

InstrumentalVariables.LinearIVModel{InstrumentalVariables.DenseIVPredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}}:

Coefficients:
     Estimate Std.Error  t value Pr(>|t|)
x1    9.89496   1.05856  9.34756   <1e-99
x2   0.280405  0.238565  1.17538   0.0836
x3   -1.27742  0.263199 -4.85346   <1e-99

To get the HC1() standard errors --- the default in STATA when ,rob is used

coeftable(iv, HC1())