How can I execute an out of sample on step forecast in the forecast package using and xreg argument in the stlm function? I am trying to program time series cross validation with the stlm function. The stlm function has the modelfunction and model arguments in order to hold the modelling constraints applied after stl decomposition constant. This should allow for me to test my model on multiple test sets. I cannot however get it to work with modelfunction = Arima whe I am also using an xreg argument. 
Below is a reproducible example. What am I doing wrong?
library(fpp)#contains sample data
library(forecast)
#Training data 
milkUse <- ts(milk[1:100], start = start(milk), frequency =  frequency(milk))
ausbeerUse <- ts(ausbeer[1:100], start = start(milk), frequency =  frequency(milk))

Model <- stlm(milkUse, "periodic", method = "arima", xreg = ausbeerUse)#inital model

#Train with different training data 
milkUse2 <- window(milkUse, start = start(milkUse), end = c(8,4))
ausbeerUse2 <- window(ausbeerUse, start = start(ausbeerUse), end = c(8,4))

stlm(milkUse2, "periodic", xreg = ausbeerUse2, modelfunction = Arima, model = Model$model) # Does not work

 A: The problem here is that stlm already has an xreg argument, so stlm is not passing it through to Arima. One solution is to redefine Arima like this:
myArima <- function(x, ...) {Arima(x, xreg=ausbeerUse2, ...)}
stlm(milkUse2, s.window="periodic", modelfunction=myArima, model=Model$model) 

So xreg is hard-coded into myArima and you don't need to pass it as an argument to stlm.
A: From gaining understanding from Rob Hyndman's answer above I have raided the forecast package github and made the below version of the stlm function. I simply removed the "auto modelling" parts of the stlm function and retained the modelfunction functionality only.
In case it helps any one else here it is.
stlm2 <- function(y ,s.window=7, robust=FALSE,
                 modelfunction=Arima, lambda=NULL, x=y, ...)
{
  # Transform data if necessary
  origx <- x
  if (!is.null(lambda))
    x <- BoxCox(x, lambda)

  # Do STL decomposition
  stld <- stl(x,s.window=s.window,robust=robust)

  # De-seasonalize
  x.sa <- seasadj(stld)

  # Model seasonally adjusted data
  fit <- modelfunction(x.sa, ...)
  fit$x <- x.sa

  # Fitted values and residuals
  fits <- fitted(fit) + stld$time.series[,"seasonal"]
  res <- residuals(fit)

  return(structure(list(stl=stld,model=fit, lambda=lambda, x=origx, m=frequency(origx),
                        fitted=fits, residuals=res),class="stlm"))
}
stlm2(milkUse2, xreg = ausbeerUse2, s.window="periodic", modelfunction=Arima, model=Model$model) 

