Combinef in R HTS package- constrain to keep forecasts positive? When using the combinef function from Rob Hyndman's very useful hts package for forecasting hierarchical and grouped time series, there does not seem to be a way to constrain the optimally combined forecasts to be positive- the starting forecasts can be positive, but can go negative through the reconciliation process.
The forecast.gts and forecast.hts functions have an argument to keep forecasts positive, but this does not seem to be an option when using combinef by itself with forecasts obtained by other methods.
Am I correct in this understanding, and if so is there a decent workaround? 
 A: The positive=TRUE argument for forecast.gts and forecast.hts ensures the starting forecasts to be positive, but not the final reconciled forecasts. Even when the starting forecasts are positive, it is possible for the reconciled forecasts to be negative. When you use combinef, you provide your own starting forecasts, so it is up to you to make them positive.
It would be possible to use a non-linear least squares reconciliation procedure to produce positively constrained reconciled forecasts, but that would be much slower.
A: You can't ensure positivity and sum consistency of hierarchical forecasts if you use forecast::combinef().
What I personally found useful was to set up the summation matrix (see the original publication by Hyndman et al., 2011) and then solve the relevant (weighted) least squares problem with additional nonnegativity constraints. This will give you nonnegative and sum consistent forecasts. I have found repeatedly that this approach still results in better forecasts on all levels in the hierarchy.
This approach also allows including equality constraints, or constraints that are more general than just "$\geq 0$". I have had applications where some forecasts needed to be larger than a certain number (because of existing orders), or needed to be constrained to be equal to a given value, which you can model using two inequality constraints.
One possible tool that solves (weighted) least squares with linear constraints is the pcls() function in the mgcv package for R. (Note that this is slightly more general than your use case: mgcv::pcls() allows for linear constraints, but you and the use cases I outline in the previous paragraph only need box constraints.) However, this is of course also not optimized to leverage the specific structure of potential forecast hierarchy matrices, so your performance may be significantly worse than if you use combinef().
A: Some code based on Rob's workaround of setting the negatives to zero and reconcile
# Re-reconciliate when zero values present

# Extract groups
groups <- hts.obj %>%
  aggts() %>%
  get_groups()

x=0 #Counter

# Loop until all positive
while(sum(hts.obj[[1]] <0) > 0){

  # Generate all time series
  hts.obj <- aggts(hts.obj)    

  # Overwrite negatives by zero
  hts.obj[hts.obj<0] <- 0

  # Reconcile
  hts.obj <- hts.obj %>%
    ts() %>%
    combinef(groups = groups, keep ="gts")

  # Count up
  x=x+1

  # Break after 10 loops
  if(x>=10)break

  }

rm("x")

# Overwrite remaining negatives by zero
hts.obj[[1]][hts.obj[[1]]<0] <- 0

In this example forecasts are indexed by [[1]] - this might change.
Also note that overwriting zeros leads to biased forecasts.
