diff --git a/NEWS.md b/NEWS.md
index 6f73733..85059e7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,9 +5,12 @@ subtitle: plm - Linear Models for Panel Data - A set of estimators and tests for
 ---
 
 # 2.6-9999 development version, changes since 2.6-4
-* vcovXX: FD models with only one observation per group prior to 
-          first-differencing errored ([#58](https://github.com/ycroissant/plm/issues/58)).
-* is.pdata.frame (non-exported helper function): fix part of the check if object 
+
+### Fixes:
+* `vcovXX`: FD models with only one observation per group prior to 
+           first-differencing errored ([#58](https://github.com/ycroissant/plm/issues/58)).
+* `pggls`: FD models errored with the data constellation as described for `vcovXX`. 
+* `is.pdata.frame` (non-exported helper function): fix part of the check if object 
                                                  does not have an index.
 
 ***
diff --git a/R/est_ggls.R b/R/est_ggls.R
index c36e18d..d450220 100755
--- a/R/est_ggls.R
+++ b/R/est_ggls.R
@@ -122,9 +122,16 @@ pggls <- function(formula, data, subset, na.action,
     if (model.name == "fd") {
     ## eliminate first year in indices
         nt <- pdim$Tint$nt[-1L]
+        n <- pdim$nT$n
         Ti <- pdim$Tint$Ti - 1
+        if(any(drop <- Ti == 0L)) {
+          # drop groups in Ti that are now empty (group had 1 observation before first-differencing, hence 0 after)
+          # and adjust n and id.names due to same reason
+          Ti <- Ti[!drop]
+          n <- n - sum(drop)
+          id.names <- id.names[!drop]
+        }
         T <- pdim$nT$T - 1
-        n <- pdim$nT$n
         N <- pdim$nT$N - pdim$Tint$nt[1L]
         time.names <- pdim$panel.names$time.names[-1L]
         groupi <- as.numeric(index[[1L]])
@@ -133,7 +140,7 @@ pggls <- function(formula, data, subset, na.action,
         sel[1L] <- 1 # the first must always be 1
         ## eliminate first obs in time for each group
         index <- index[!sel, ]
-        id <- index[[1L]]
+        id <- droplevels(index[[1L]])
         time <- factor(index[[2L]], levels = attr(index[ , 2L], "levels")[-1L])
     } else {
         nt <- pdim$Tint$nt
@@ -229,7 +236,7 @@ pggls <- function(formula, data, subset, na.action,
         }
 
         subOmega <- rowMeans(tres, dims = 2L, na.rm = TRUE) # == apply(tres, 1:2, mean, na.rm = TRUE) but faster
-        list.cov.blocks <- sapply(other.list, function(i) subOmega[i, i], USE.NAMES = FALSE)
+        list.cov.blocks <- sapply(other.list, function(i) subOmega[i, i], USE.NAMES = FALSE, simplify = FALSE)
         omega <- bdsmatrix::bdsmatrix(groupsdim, unlist(list.cov.blocks, use.names = FALSE))
     }
     A <- crossprod(X, solve(omega, X))
diff --git a/inst/tests/test_FD_models.R b/inst/tests/test_FD_models.R
index c92e2e1..1ffbcb2 100644
--- a/inst/tests/test_FD_models.R
+++ b/inst/tests/test_FD_models.R
@@ -49,9 +49,10 @@ mod <- plm(inv ~ value + capital, data = pGrun1, model="fd")
 vcovHC(mod)
 vcovBK(mod)
 diff(pGrun1$inv)
+ 
+## pggls run test
+mod.pggls <- pggls(inv ~ value + capital, data = pGrun1, model="fd")
 
-## TODO: pggls errors
-# mod.pggls <- pggls(inv ~ value + capital, data = pGrun1, model="fd")
 
 # one time period with only one observation
 pGrun2 <- pGrun[-c(21, 41, 61, 81, 101, 121, 141, 161, 181), ]
diff --git a/inst/tests/test_FD_models.Rout.save b/inst/tests/test_FD_models.Rout.save
index 09deaf8..a7419cf 100644
--- a/inst/tests/test_FD_models.Rout.save
+++ b/inst/tests/test_FD_models.Rout.save
@@ -322,9 +322,10 @@ attr(,"cluster")
   11.9   32.2  -32.6    3.5   26.3   38.6  -21.1  -30.6   -4.5   36.8   66.3 
 3-1947 3-1948 3-1949 3-1950 3-1951 3-1952 3-1953 3-1954 
  -12.7   -0.9  -48.0   -4.8   41.7   22.1   22.2   10.1 
+>  
+> ## pggls run test
+> mod.pggls <- pggls(inv ~ value + capital, data = pGrun1, model="fd")
 > 
-> ## TODO: pggls errors
-> # mod.pggls <- pggls(inv ~ value + capital, data = pGrun1, model="fd")
 > 
 > # one time period with only one observation
 > pGrun2 <- pGrun[-c(21, 41, 61, 81, 101, 121, 141, 161, 181), ]
@@ -363,4 +364,4 @@ attr(,"cluster")
 > 
 > proc.time()
    user  system elapsed 
-   1.59    0.20    1.78 
+   1.90    0.23    2.17