How to estimate a fixed effects regression WITH robust standard errors AND instrument variables I have been trying to estimated the stated problem, but I only succeed in parts of it. The following artificial setup is supposed to illustrate my problem in detail:
Setup the data:
library(plm)
data(Cigar)

Fixed effects AND robust standard errors  (works perfectly)
  fit <- plm(price ~ sales + pop, data=Cigar, index=c("state","year"), model="within")

> coeftest(fit, vcov.=vcovHC(fit))

t test of coefficients:

        Estimate Std. Error t value  Pr(>|t|)    
sales -1.0391141  0.1671141 -6.2180 6.726e-10 ***
pop    0.0190151  0.0064447  2.9505  0.003228 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

New Fixed effects regression WITH Instruments: (sales is instrumented with cpi)...
fit <- plm(price ~ sales + pop | cpi + pop, data=Cigar, index=c("state","year"), model="within")

...renders usage of coeftest problematic!

coeftest(fit, vcov.=vcovHC(fit))
      Fehler in vcovG.plm(x, type = type, cluster = cluster, l = 0, inner = inner,  : 
        Method not available for IV

Well the error message is obvious, but I am just wondering what other people do who face the same problem as I do?
I hope I don't miss something but there is no statistical reason why heteroskedasticity cannot apply in this setting, right?
So is there any possibility to estimate robust standard errors in this setting with the plm package or are there any R-packages which may help? 
Thank you very much for your answers - any advice is appreciated!
 A: EDIT: The methods are also in CRAN versions of plm >= 1.6-4.
The appropriate methods for robust vcovs were not implemented in plm. They are now in the development version >= 1.6-1 (see http://r-forge.r-project.org/R/?group_id=406).
library(plm)
library(lmtest)
data(Cigar)
fit <- plm(price ~ sales + pop, data=Cigar, index=c("state","year"), model="within")
fit2 <- plm(price ~ sales + pop | cpi + pop, data=Cigar, index=c("state","year"), model="within")

coeftest(fit, vcov.=vcovHC(fit))

# t test of coefficients:
#
#
#         Estimate Std. Error t value  Pr(>|t|)    
# sales -1.0391141  0.1671141 -6.2180 6.726e-10 ***
# pop    0.0190151  0.0064447  2.9505  0.003228 ** 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

coeftest(fit2, vcov.=vcovHC(fit2))

# t test of coefficients:
#
#        Estimate Std. Error t value Pr(>|t|)   
# sales -6.2479556  1.9032780 -3.2827 0.001055 **
# pop   -0.0021752  0.0121773 -0.1786 0.858260   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

