Optimal combination approach and forecasted proportions

Question

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?)

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.