diff --git a/ChangeLog b/ChangeLog
index c0cdf23..421b368 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,7 +1,8 @@
-v4.4 (xx 2014)
+v4.4 (17 July 2014)
 * Allowed multiple seasonal objects "msts" in hts() and gts().
 * Fixed bug of MASE in accuracy.gts().
 * Allowed user's defined forecasting function in forecast.gts().
+* Fixed bug of xreg and newxreg in forecast.gts(), when a matrix is passed.
 
 v4.3 (10 June 2014)
 * Fixed bug of the arg "include" in plot.gts, when there are not yearly data.
diff --git a/DESCRIPTION b/DESCRIPTION
index 0eea6ba..ff99152 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: hts
 Type: Package
 Title: Hierarchical and grouped time series
-Version: 4.3
+Version: 4.4
 Depends: forecast (>= 5.0)
 Imports: SparseM, parallel, utils
 Suggests: testthat
diff --git a/R/forecast.gts.R b/R/forecast.gts.R
index 6b378a1..b0169fa 100644
--- a/R/forecast.gts.R
+++ b/R/forecast.gts.R
@@ -112,24 +112,66 @@ forecast.gts <- function(object, h = ifelse(frequency(object) > 1L,
     return(out)
   }
 
-  if (parallel) { # parallel == TRUE
-    if (is.null(num.cores)) {
-      num.cores <- detectCores()
+  if (fmethod == "arima" && !is.null(c(xreg, newxreg)) 
+      && !is.vector(xreg) && !is.vector(newxreg)) {
+    # Error handling for args xreg and newxreg
+    nc.y <- ncol(y)
+    nr.y <- nrow(y)
+    # the dim of xreg and newreg should be equal to the dim of forecasting ts 
+    check.col <- all(nc.y == c(ncol(xreg), ncol(newxreg)))
+    check.row <- nr.y == nrow(xreg)
+    check.row.new <- h == nrow(newxreg)
+    if (!check.col) {
+      msg <- sprintf("The args xreg or newxreg should have %i columns.", nc.y)
+      stop(msg)
+    } else if (!check.row) {
+      msg <- sprintf("The arg xreg should have %i rows", nr.y)
+      stop(msg)
+    } else if (!check.row.new) {
+      msg <- sprintf("The arg newxreg should have %i rows", h)
+      stop(msg)
+    }
+    # parallel process not working and conventional forloop will be used instead
+    # because auto.arima and forecast.Arima use the same arg xreg, and lapply
+    # cannot be able to tell xreg and newxreg
+    parallel <- FALSE 
+    pfcasts <- matrix(, nrow = h, ncol = nc.y)
+    if (keep.fitted) {
+      fits <- matrix(, nrow = nr.y, ncol = nc.y)
+    }
+    if (keep.resid) {
+      resid <- matrix(, nrow = nr.y, ncol = nc.y)
+    }
+    for (i in 1L:nc.y) {
+      models <- auto.arima(y[, i], lambda = lambda, xreg = xreg[, i], ...)
+      pfcasts[, i] <- forecast(models, h = h, xreg = newxreg[, i])$mean
+      if (keep.fitted) {
+        fits[, i] <- fitted(models)
+      }
+      if (keep.resid) {
+        resid[, i] <- residuals(models)
+      }
+    }
+  } else {  # if xreg and newxreg are not used
+    if (parallel) { # parallel == TRUE
+      if (is.null(num.cores)) {
+        num.cores <- detectCores()
+      }
+      cl <- makeCluster(num.cores)
+      loopout <- parSapplyLB(cl = cl, X = y, FUN = function(x) loopfn(x, ...), 
+                             simplify = FALSE)
+      stopCluster(cl = cl)
+    } else {  # parallel = FALSE
+      loopout <- lapply(y, function(x) loopfn(x, ...))
     }
-    cl <- makeCluster(num.cores)
-    loopout <- parSapplyLB(cl = cl, X = y, FUN = function(x) loopfn(x, ...), 
-                           simplify = FALSE)
-    stopCluster(cl = cl)
-  } else {  # parallel = FALSE
-    loopout <- lapply(y, function(x) loopfn(x, ...))
-  }
 
-  pfcasts <- sapply(loopout, function(x) x$pfcasts)
-  if (keep.fitted) {
-    fits <- sapply(loopout, function(x) x$fitted)
-  }
-  if (keep.resid) {
-    resid <- sapply(loopout, function(x) x$resid)
+    pfcasts <- sapply(loopout, function(x) x$pfcasts)
+    if (keep.fitted) {
+      fits <- sapply(loopout, function(x) x$fitted)
+    }
+    if (keep.resid) {
+      resid <- sapply(loopout, function(x) x$resid)
+    }
   }
 
   if (is.vector(pfcasts)) {  # if h = 1, sapply returns a vector
diff --git a/man/forecast.gts.Rd b/man/forecast.gts.Rd
index e514b2f..3c2c576 100644
--- a/man/forecast.gts.Rd
+++ b/man/forecast.gts.Rd
@@ -36,8 +36,8 @@ Methods for forecasting hierarchical or grouped time series.
   \item{FUN}{It allows users to define their own fitted model that can be passed to the \code{forecast} function
     to generate base forecasts such as \code{tbats} and \code{stlf}. There is no use in the 
     arguments \code{positive} and \code{lambda}, if \code{FUN} is being used.}
-  \item{xreg}{When \code{fmethod = "arima"}, a vector or matrix of external regressors, which must have the same number of rows as the original univariate time series}
-  \item{newxreg}{When \code{fmethod = "arima"}, a vector or matrix of external regressors, which must have the same number of rows as the original univariate time series}
+  \item{xreg}{When \code{fmethod = "arima"}, a vector or matrix of external regressors used for modelling, which must have the same number of rows as the original univariate time series}
+  \item{newxreg}{When \code{fmethod = "arima"}, a vector or matrix of external regressors used for forecasting, which must have the same number of rows as the \code{h} forecast horizon}
   \item{...}{Other arguments passing to \code{\link[forecast]{ets}} or \code{\link[forecast]{auto.arima}}}
 }
 \value{A forecasted hierarchical/grouped time series of class \code{gts}.}