I referred to this link and I have the following questions regarding my data. Let me start by explaining the time series that I am dealing with.

I have daily hospital data with various departments and numerous doctors working in each department. I have several years of data and my forecast horizon is for the next 365 days. My data has weekly and annual seasonality. Moreover I intend to capture the effects of holidays and Sundays in my forecasts. As a result I have not created a hierarchical time series as suggested towards the end of the link(primarily because I am not sure whether we can pass a regressor to it and more so because I do not know how many doctors I end up predicting for in each department).

The reason for this is that some doctors do not have good data(short time series or sparse data). In this case I collect these doctors and aggregate them to form something I call "OtherDocs". Typically in DeptXYZ -> Doc1 , Doc2 , Doc3 , Doc4 , Doc5 and Doc6 I could end up creating forecasts for DeptXYZ -> Doc1 , Doc3 , Doc4 , Doc6 and OtherDocs. If OtherDocs is still not predictable I generate a naive forecast. In this fashion I created base forecasts for every level in the hierarchy individually using arima and passing my xreg to it and selecting the best model on the basis of AIC.

Now, consider this example -

Total -> DeptX and DeptY

DeptX -> DocA and DocB

DeptY -> Doc1 , Doc2 and Doc3

There are cases where DocA has a time series that starts from "2011-03-11" and ends on "2016-09-07" while DocB has a time series that starts from "2011-05-17" and ends on "2016-09-07". Generating the base forecasts for DocA and DocB results in the predicted values(fit$mean) being of a time series from "2016-09-08" to "2017-09-07". As long as the time series refers to the same dates within the Department I believe we are good to go.

In my attempt to reconcile the forecasts from each level I employed the forecasted proportions like so -

$\Largeỹ_{DocA,365} = \frac{ŷ_{DocA,365}*ŷ_{DeptX,365}}{(ŷ_{DocA,365}+ŷ_{DocB,365})*(ŷ_{DeptX,365}+ŷ_{DeptY,365})}ŷ_{Total,365}$

1. Am I doing anything wrong in the above step?

2. Suppose for one moment that the topmost level forecasted values do not capture the low points of data in the case of Holidays and Sundays. Does that intuitively mean that revised forecasts for DocA might not correctly capture the same(being a proportion of $ŷ_{Total,365}$)?

Another query I have is to do with the Optimal Combination Approach -

$\Largeỹ_h = S(S′S)^{-1}S′ŷ_h$

3. I am unfamiliar with this matrix notation $S'$. Is it the inverse of $S$? Could you shed some light on this? And how do you suggest I calculate the summing up matrix in my case?(Is it absolutely necessary to proceed with the exact knowledge of the number of doctors in each department?)

  • $\begingroup$ I think this $S'$ implies this $S^T$. As in, it is the transpose of the matrix. $\endgroup$
    – Ic3fr0g
    Sep 10, 2016 at 8:05

1 Answer 1


In response to your 3 questions:

  1. Yes, that notation looks correct to me.
  2. In general, no, if your top level time series does not capture a feature of your data that manifests at the top level, your lower level series won't either. It is possible that if one of your doctors capture the dip but the rest don't, that doctor will see a reduced forecast (since it's forecast is a smaller relative to the other doctors), but since all the values need to add up to the top level forecast, you may actually see an increase in the forecast for the rest of the doctors in response. In general the reason to do a top down style of forecast is that higher level forecasts tend to smooth out randomness and do a better job capturing complexities in your data. If this is not the case in your data, top down may not be for you.
  3. Mayur H is correct that S' is indicating a transpose operation. As for calculating the summing matrix, by the time you are aggregating forecasts, you should know how many forecasts you have to aggregate. This is something you need to know before aggregation, but I am unclear what the circumstance might be that would prevent you from having this information by the time you had decided which series you are forecasting.

As for your decision not to use the R hts package, you may want to reconsider. You are still able to make ARIMA forecasts with external regressors by making the function call:

forecast(htsObject, h = horizon, method = tdfpORcomb, fmethod = "arima",
         xreg = trainingRegressors, newxreg = forecastingRegressors)

The first two arguments should be self explanatory, method would be "tdfp" for top down forcasted proportions or "comb" for the optimal combination, xreg are the external regressors for your training data and newxreg would the those for your forecast time points. As for your concerns about the poor quality of data for some doctors, the hts package can handle that. Even very short time series can be forecast using this method since information from other parts of the hierarchy can be used to inform short time series forecasts. Lets take as an example below. For our time series, we will do the following

> set.seed(1)
> nodes <- list(2, c(3, 2))
> abc <- ts(5 + matrix(sort(rnorm(100)), ncol = 5, nrow = 20), 
            frequency = 4)
> abc[1:18, 3] <- NA
> x <- hts(abc, nodes)
> x$bts
     Series 1 Series 2 Series 3 Series 4 Series 5
1 Q1 2.785300 4.387974       NA 5.387672 5.782136
1 Q2 3.010648 4.410479       NA 5.389843 5.821221
1 Q3 3.195041 4.426735       NA 5.398106 5.881108
1 Q4 3.476433 4.431331       NA 5.417942 5.918977
2 Q1 3.529248 4.457480       NA 5.475510 5.943836
2 Q2 3.622940 4.521850       NA 5.487429 6.063100
2 Q3 3.723408 4.526599       NA 5.556663 6.100025
2 Q4 3.746367 4.556708       NA 5.558486 6.124931
3 Q1 3.775387 4.585005       NA 5.569720 6.160403
3 Q2 3.870637 4.605710       NA 5.575781 6.178087
3 Q3 3.955865 4.632779       NA 5.593901 6.207868
3 Q4 4.065902 4.694612       NA 5.593946 6.358680
4 Q1 4.164371 4.695816       NA 5.610726 6.433024
4 Q2 4.179532 4.746638       NA 5.619826 6.465555
4 Q3 4.256727 4.835476       NA 5.689739 6.511781
4 Q4 4.290054 4.844204       NA 5.696963 6.586833
5 Q1 4.292505 4.864821       NA 5.700214 6.595281
5 Q2 4.311244 4.864945       NA 5.738325 6.980400
5 Q3 4.373546 4.887654 5.364582 5.763176 7.172612
5 Q4 4.378759 4.897212 5.370019 5.768533 7.401618

We now have 5 time series, but one of them (Series 3) is too short to forecast using ARIMA.

> auto.arima(x$bts[ , 3])
Error in OCSBtest(x, m) : The OCSB regression model cannot be estimated

However, because of the additional information in the hierarchy, we can get a forecast for Series 3 using the hts version of the forecast function:

> xForecast <- forecast(x, h = 10, method = "comb", fmethod = "arima")
> xForecast$bts
     Series 1 Series 2  Series 3 Series 4 Series 5
6 Q1 4.406996 4.922615  5.295991 5.770775 7.573794
6 Q2 4.430893 4.951091  5.384539 5.798779 7.632580
6 Q3 4.458601 4.978792 10.796670 5.813754 7.750316
6 Q4 4.482059 5.005933 10.820060 5.834600 7.836871
7 Q1 4.503295 5.031758 10.768352 5.849883 7.941606
7 Q2 4.530396 5.060007 10.844907 5.870614 8.038023
7 Q3 4.556540 5.087792 16.261433 5.888068 8.142945
7 Q4 4.580688 5.114878 16.281914 5.907109 8.238842
8 Q1 4.601632 5.140740 16.232187 5.923707 8.336769
8 Q2 4.628847 5.168963 16.307364 5.943480 8.438147

Before investing too much time rewriting the logic of the hts package, it may be worth trying it to see if you can save time and get the results you want.

  • $\begingroup$ I did end up rewriting the TDFP and Combination Approach. Thanks for a great answer. I will implement this soon and get back to you. I realised that I couldn't capture the dips during holidays because the revised_doc_forecasts are a proportion of the top_level aggregated forecasts. I intend to apply the middle out approach in my case. For some reason I can't upvote your answer(need 15 rep). $\endgroup$
    – Ic3fr0g
    Sep 14, 2016 at 5:03
  • 1
    $\begingroup$ Please upvote my question so that I can in turn upvote your answer! $\endgroup$
    – Ic3fr0g
    Sep 14, 2016 at 6:12
  • $\begingroup$ If you are having issues capturing the effects of holidays and Sundays, there are two things that may be useful. First, give your series a frequency of 7 to allow the model to find the "Sunday" effect. Second, add a "holiday" external regressor that is 1 for holidays and 0 otherwise. As for your method, the best choice depends on your data. Middle out may help you if your "Dept" level is dead on, and the rest are worse, but most of the time it isn't that clear cut. Try using comb with weights = "sd", as this will take into account model fit when redistributing the forecasts. $\endgroup$
    – Barker
    Sep 14, 2016 at 17:08
  • $\begingroup$ Ofcourse the forecasted data at the dept level looks solid. But not at the aggregated level. But I realised that I was going wrong with my implementation of Optimal Combination. It has been fixed and solved. Thanks once again $\endgroup$
    – Ic3fr0g
    Sep 15, 2016 at 5:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.